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ABSTRACT 

Observed spectral energy distributions (SEDs) of FU Orionis, V1057 Cygni, and 
V1515 Cygni are fit by theoretical spectra, which are calculated from models consisting 
of outbursting accretion disks together with flattened envelopes. Temperature in 
the envelopes is determined by approximate radiative equilibrium with a central 
source. The disk models are two-dimensional and include reprocessing of disk 
radiation by the disk. The theoretical spectra are calculated using a radiative transfer 
code and frequency-dependent opacities, at a spectral resolution of A/AA = 14. 
Excellent matches to the data are obtained for all three objects with reasonable model 
parameters. Radiative transfer is also used to calculate a time series of images showing 
the progress of an outburst as imaged through a B-band filter. 


Subject headings: accretion disks: protostellar — stars: pre-main-sequence — stars: 
FU Orionis — stars: V1057 Cygni — stars: V1515 Cygni 


1. Introduction 

The FU Orionis phenomenon (reviewed by Herbig 1966, 1977; Hartmann, Kenyon & Hartigan 1993) 
represents a key event in the early history of a star which should provide us with important clues 
on the nature of pre-main-sequence stellar evolution and the formation of planetary systems. The 
three objects considered here are FU Orionis, which brightened by 6 magnitudes in 120 days in 
1936-7, V1057 Cygni, which brightened by 5.5 magnitudes in 250 days in 1969, and V1515 Cygni, 
which showed a much longer rise time of about 20 years, reaching a peak about 1970. At present 
the three objects are respectively 1 mag, 3 mag, and 0.5 mag fainter than at maximum light. A 
recent observational summary is provided by Bell et al. (1995), henceforth referred to as BLHK. 


^CO/Lick Observatory Bulletin No. ... 



- 2 - 


The most likely explanation for the observed outbursts is rapid mass accretion onto the surface 
of a star, triggered by a thermal instability in the inner, ionized regions of a surrounding accretion 
disk (Paczynski 1976; Hartmann & Kenyon 1985; Lin & Papaloizou 1985). The principal physical 
effect that induces the instability is a heating rate, say by viscous dissipation, that increases more 
rapidly with temperature than does the cooling rate by radiative and convective transfer. This 
situation is likely to occur in the region of partial ionization of hydrogen where the opacity is a 
strongly increasing function of temperature. The details of the mechanism remain to be clarified, 
since current models rely on an arbitrary viscous dissipation parameter (a) and in some cases 
an arbitrary external perturbation. Objections have been raised, on observational grounds, to 
the disk model (Herbig 1989; Petrov & Herbig 1992), but by now it is clear that steady-state 
disk models can explain a number of the observed features (Hartmann & Kenyon 1985; Adams, 
Lada, & Shu 1987; Kenyon, Hartmann, & Hewett 1988; Kenyon & Hartmann 1989, 1991; Calvet, 
Hartmann, & Kenyon 1993; Hartmann, Kenyon, & Hartigan 1993). Moreover, certain features 
which are still difficult to explain with a model involving a disk with constant mass flux can be 
explained by the use of more elaborate, time-dependent, vertically resolved models (Kawazoe & 
Mineshige 1993; Bell & Lin 1994; BLHK). 

The detailed models of BLHK, which take into account the radial and vertical structure of 
the disk as a function of time through a series of outbursts, form the basis of the present paper. 
The general picture is that the star-disk system is still at a sufficiently early stage of its evolution 
so that matter is falling onto the outer part of the disk from the surrounding interstellar cloud 
(Kenyon & Hartmann 1991) at a rate M{ n . The results of the calculations show that as long as 
M( n exceeds 5xl0 -7 M 0 yr” 1 , independent of the viscosity parameter a, the thermal ionization 
instability is initiated only a few stellar radii out from the star, and the disk goes into outburst. 
The ionization front then propagates radially through the disk on a time scale a few times 
longer than the local thermal value, reaching a maximum radius not larger than 0.25 AU (Bell & 
Lin 1994). The duration of the outburst is determined by the local viscous diffusion time in the 
hot, ionized region. Once this region has been sufficiently depleted by accretion onto the central 
star, it cools, hydrogen recombines, and the outburst declines. The interval between outbursts is 
determined by M* n , which controls the rate at which mass is added to the inner disk. Once the 
mass there has increased to the point where the temperature once again becomes high enough to 
ionize hydrogen, the outburst cycle repeats. Spontaneous outbursts by this mechanism can explain 
the light curve of V1515 Cyg, but to explain the more rapid rise times of FU Ori and V1057 Cyg 
a small density perturbation in the inner disk is required to initiate the instability. This general 
explanation has been used by Lin et ai (1994) to conclude that in the case of HL Tau the value of 
M{ n is much larger than the accretion luminosity implied from the properties of the inner disk, so 
the object could be in the quiescent stage between outbursts. 

The work of BLHK emphasizes comparison of models with observed light curves. In this 
paper we make further tests of the viability of the detailed two-dimensional models by comparing 
them with observed spectral energy distributions and photometric properties. Previous work 
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on comparison of disk models of the overall spectral energy distribution (Adams et al. 1987, 
Kenyon et al. 1988; Kenyon k Hartmann 1991; BLHK) have assumed that at each radius the 
theoretical disk radiates either the same spectrum as a black body at the disk surface temperature, 
or as a stellar photosphere at the same effective temperature. A different approach is followed 
here. Frequency-dependent radiative transfer models are obtained from the two-dimensional disk 
structures (§ 2). Reprocessing of light from the hot inner disk by the cooler outer disk surface 
is included. Theoretical spectra are obtained as a function of the angle of inclination between 
the rotation axis and the line of sight (§3) and the differences between spectra that include 
reprocessing and those that don’t are emphasized. The results are compared in detail (§ 4) with 
observations of the spectral energy distributions of the three objects in the wavelength range 
0.3—100 pm. The frequency-dependent information is then used to provide a detailed calculation 
of the light distribution as a function of position and time for the model that fits the light curve 
of V1515 Cyg. The implications of the results are discussed in § 5. 


2. Numerical Disk models 

In this section, we describe the disk models used to investigate the appearance of FU Orioms 
systems. The first component, outlined in § 2.1, is a one-dimensional, radial time-dependent set 
of disk evolution equations with parameters chosen in BLHK to produce an outburst with B light 
curves close to those observed in one of the three program objects FU Orionis, V1515 Cygni, 
and V1057 Cygni. A particular epoch of interest is then selected from the solution to the 
time-dependent equations, and the model is expanded (§ 2.2) into a set of vertical structure models 
at a sufficient number of radii to define T{r,z) and p{r,z ), where r is the cylindrical radius and 
z is height above the midplane. Reprocessing caused by the illumination of one part of the disk 
surface by another is found to significantly raise surface temperatures at some radii. Our method 
for calculating this reprocessing is described in § 2.3. When the disk is quiescent and its luminosity 
is low, the central star may make an important contribution (§ 2.4). Finally in § 2.5, we add an 
envelope to account for the observed flat spectrum at long wavelengths. The resulting model disk, 
consisting of the vertical models modified by reprocessing, plus an envelope and central star, is 
shown schematically in Figure 1. From this structure, the spectra presented in § 4 are calculated 
as discussed in § 3. 


2.1. Time Dependent Models 

The outbursting disk is evolved through time by a numerical integration of a set of equations 
which represent the radial diffusion of disk material by viscosity. The equations are described in 
detail by Bell k Lin (1994) and consist of mass continuity, energy balance, and the (f> component 
of the equation of motion (Bell k Lin, eqns 3-5). The energy equation includes viscous heating 
from a Shakura k Sunyaev (1973) a viscosity, radiative losses from the surface, radial diffusion of 
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Fig. L— Model FU Orionis objects consist of a disk D and an envelope E. The envelope has a 
thickness A z and a central hole of radius R^, The object is viewed at an angle i from pole-on. The 
central star is located at the lower left. Also shown are three example lines-of-sight along which 
radiative transfer is carried out. 




-5- 


radiation, and PdV work. In the region of outburst, the disk is not necessarily in vertical thermal 
balance. The opacity used in this initial stage is the analytic approximation to the Rosseland 
mean opacity described in the Appendix of Bell & Lin (1994). The values of a are the same as 
those used by BLHK: 10“ 4 throughout most of the disk where hydrogen is neutral, and 10“ 3 in 
regions where the hydrogen is ionized. 

BLHK found that the observed light curves of the three objects were fit acceptably over 
a range of model parameters; these parameters include Mi n , and the size and type of density 
perturbation to the inner disk. In this paper we choose their models Al, Bl, and Cl from which 
to calculate detailed spectral energy distributions; the models’ parameters are listed in Table 1. 
Outburst in the Bl disk is spontaneous, but in Al and Cl a small perturbation is required in order 
to match the rapid rise times of the FU Ori and V1057 Cyg outbursts. Such a perturbation might 
be provided in a clustered star-forming environment by interactions with nearby protostars. All 
three program objects lie in associations of young stars (see §3.4). Alternatively, the perturbations 
might be caused by radial migration of protoplanets (Syer & Clarke 1996) or by time variation of 
the effective viscosity in the outer disk. The perturbations are annuli of material added to the 
radial time-dependent models during the quiescent phase between outbursts. They have inner 
radii r p , masses M p , and a fractional surface density change within the annulus of AS/E. Outer 
radii are determined by M p and AS/E. 


2.2. Vertical Structure Models 

In this subsection we discuss the sets of vertical structure models used to calculate spectra 
and images of the disk. The models provide the detailed distributions of temperature and density 
throughout the disk, T(r,z ) and p(r,z). The vertical structure procedure is described in detail in 
Bell & Lin (1994); the essence of the program and several modifications are described here and 
in § 2.3. The basic procedure is as follows: 

1. Start by calculating a set of vertical structure models, for example, for model Cl of BLHK. 


Table 1: Disk outburst model parameters 


Star 

Model 

Min 

M 0 /yr 

M P / 

Mq 

AS/S 

r vl 

r Q 

FU Orionis 

Al 

3 x 1CT 6 

0.01 

5 

13 

V1515 Cygni 

Bl 

1 x 1(T 5 

0 

0 

— 

V1057 Cygni 

Cl 

1 x 10" 6 

0.002 

3 

10 


Note: Table 2 in Bell et al 1995 (BLHK) contains two misprints, which are corrected here. 
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2. Use the resulting disk surface properties, T e (r) and H(r), to calculate the flux of energy 
deposited at a given point by illumination from other points of the disk. 

3. Add this flux at the level where Rosseland optical depth f = 2/3 to get a new effective 
temperature at each radius. 

4. Calculate a new set of vertical structure models with this new boundary condition. 

In the remainder of this subsection we discuss the first of these steps. The remaining three points 
are discussed in § 2.3. 

The radial solution provides the disk effective surface temperature T e and degree of departure 
from vertical thermal balance, at each radius. Calculation of the vertical structure model begins 
near the surface of the disk. Surface pressure is chosen so as to balance the weight of the overlying 
atmosphere, which is otherwise ignored. In the vertical models calculated before reprocessing, 
the boundary lies at r = 2/3, while in the final vertical models from which spectra are computed 
(§ 2.3), it lies at r = 0.03. Pressure and flux are next integrated towards the disk midplane. The 
pressure gradient balances the exact expression for the gravitational attraction of a central point 
mass - that is, finite disk thickness H « r is allowed. The flux gradient includes terms due to local 
energy generation by viscous dissipation, and transport by both mixing-length-theory convection 
and radiative diffusion. Diffusion in the optically thin surface layer (and throughout the vertical 
model) is flux limited as in Bodenheimer et al. (1990), and scattering is assumed to be negligible. 
The opacity is the same approximation to the Rosseland mean as used by Bell & Lin (1994), 
although the final vertical models use the more complete opacities discussed in § 2.3. Finally, the 
disk height is adjusted and the integration is repeated until it yields zero flux at the midplane. 

Vertical structure models are calculated at 190 values of the radius, evenly-spaced in log r from 
r = 3 Rq to r = 158 AU. Time dependent models provide radial points out to 100R© fs 0.5 AU. 
Beyond this, the disk is assumed to be transporting mass in vertical thermal balance at the rate 
M in . Vertical models each consist of 100 points, with a grid adaptively-spaced so as to resolve 
steep z gradients in temperature occurring near the surface: the gap between points is held to no 
more than one twentieth of the local temperature scale height T/(dT/dz). 


2.3. Reprocessing 

During outbursts, models’ luminosities may exceed 100 L Q , most of which is emitted by the 
portion of the disk inside 0.25 AU (the “inner disk”). Some fraction of the photons emitted by 
the inner disk is intercepted by the outer disk, heating it significantly. While the inner disk has a 
spectrum peaking in the visible, the outer disk is heated to temperatures such that its peak flux 
is at wavelengths A > 2 /an. The additional flux proves useful in matching the observed fluxes at 
wavelengths 2-5 pm (§ 4.1). The effect of this “reprocessing” has been studied in the limit of a 
thin disk with central star by Adams & Shu (1986). In comparison, a disk whose surface is concave 
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up, or flared, occupies a larger solid angle as viewed from the central source, and intercepts a 
larger fraction of the photons (Kenyon & Hartmann 1987; Ruden & Pollack 1991). All of these 
calculations, however, approximate the central source by a sphere of uniform temperature. 

The construction of a pseudo-two-dimensional disk, described in the previous two subsections, 
results in a disk surface with several local maxima as shown in Figure 2. In calculating the 
temperature distribution arising from reprocessing in this case, one must consider how the wrinkles 
hide the bright inner disk from some regions of the outer disk. The details of the numerical 
procedure developed to handle this problem will be discussed by Bell (1997). The calculation 
assumes that the disk is solid at the t — 2/3 surface, and radiates and absorbs like a black body. 
The calculation is done for the entire disk, not accounting for the presence of an envelope. The 
disk surface is divided into concentric rings by cylindrical radius cu, and the rings are subdivided 
by equatorial angle (f> into surface elements. To calculate the reprocessing at a point (u>o,</> = 0) 
on the disk surface, a line of sight is drawn to each other point (u>, <f>). Normals to the surface at 
the two points are projected onto the line of sight between them. If the components are directed 
towards one another, there is possible mutual heating. Finally, the line of sight is traced to 
determine whether an intervening fold of the disk obscures the end points from one another. The 
sum Fi n of the fluxes arriving from all emitting points determines the reprocessing temperature 
T rp at the point (a> 0 ,0), in balance with losses by black body radiation, via 

C tT t 4 p = Fin. (!) 

The revised effective temperature at (a>o,0) is then 

T e ,tot = (T e 4 + T^) 1 / 4 . (2) 

New vertical structure models are next computed starting from the revised effective 
temperatures. Since the surface layer will be important in the radiative transfer calculation, the 
vertical structure calculation is now begun at an optical depth of to = 0.03. The temperature 
at this new starting point, T 0 , is found from the effective temperature using the Eddington 
approximation, T 0 4 = |T 4 fot (^ + 2/3). If the resulting T 0 is less than the reprocessing temperature 
T rp , it is set equal to T rp . This approximately represents the deposition of energy into the disk 
atmosphere by reprocessing (calculated in detail by Calvet et al. 1991), which was ignored in the 
calculation of lines-of-sight across the r = 2/3 surface. Because the flux to be reprocessed enters 
the disk from both above and below, the net flux through the disk surface is unchanged when 
reprocessing is included and in vertical thermal balance is still equal to the flux due to viscous 
energy generation. Thus the surface flux used in the calculation of the new vertical structure is 
still oTg. 

The Rosseland mean opacities used in the new vertical structure models were provided 
by Alexander (1995); they correspond to the frequency-dependent opacities described in § 3.1 
and used in the radiative transfer calculation of § 3.2. At temperatures and densities outside 
the table calculated by Alexander (1995), the Rosseland mean opacities of Pollack et al. (1985) 
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(low temperatures), Alexander, Auguson, & Johnson (1989) (higher temperatures), and Cox & 
Tabor (1976) (any remaining points) are used. 

Temperature profiles for the B1 model with and without reprocessing are shown in Figure 3. 
Without reprocessing, the temperature outside the outburst region varies as r~ 3 / 4 , as in standard 
constant mass-flux accretion disk models. With reprocessing, there is a slight enhancement in 
temperature in the outburst region, inside the radius where the disk reaches a local maximum in 
thickness. Here, opposite sides of the disk are tilted towards one another across the central axis. 
We refer to this region as the “volcano”. The most important effect of reprocessing, however, 
occurs at radii 0.1-10 AU, wherever the disk scale height is large enough so the surface sees over 
the near rim of the volcano to the hot region within. Because the reprocessing is set by the run 
of disk scale height with radius, and the scale height increases with opacity (Bell et ai 1997), 
the inner edge of this reprocessing region lies near the minimum in the disk thickness, which 
corresponds to the radius where the temperature is just low enough for dust grains to condense. 
In the spectra calculated in § 4.1, only the inner edge of the reprocessing region plays a role. 
The rest is hidden beneath an envelope. However, the radiative equilibrium approximation used 
in the envelope (§ 2.5) yields a temperature nearly continuous with that of the disk beneath, 
which was calculated by reprocessing for the disk alone. The match indicates that the disk surface 
temperatures are also approximately correct for the case where the disk is covered by the envelope. 
Because the envelope is optically thick, the surface temperature of the disk beneath the envelope 
plays no part in determining the SED. 

Over the range of radii where reprocessing is important, the disk’s temperature profile is 
flattened almost to T ~ r -1 / 2 . As one moves from the inner edge of this range to the outer, the 
disk height increases faster than the radius, the surface turns concave away from the midplane or 
“flared”, and an increasing fraction of the inside of the volcano is visible from the disk surface. 
Beyond r « 10 AU, reprocessing is insignificant because the disk surface is concave towards the 
midplane and the bright inner disk is hidden from view. This change in curvature is due mostly 
to the opacity law k(T) (Bell et ai 1997). Below the water ice condensation temperature of 125 K 
to which the disk material falls near 10 AU, k is proportional to T 2 , while in warmer material the 
opacity varies more slowly with temperature. 

The main change in the disk’s vertical structure due to inclusion of reprocessing, apart from 
the rise in surface temperature, is an increase in the thickness of the surface layer. The increase 
is greatest at radii between 4 and 8 AU, where reprocessing thickens the disk by 10%. At the 
midplane, temperature and density nowhere change more than 1.6%. For this reason, the extra 
energy input from reprocessing is not expected to significantly change the time evolution of the 
radial disk model. The calculation of reprocessing and resulting new vertical structure models 
can in principle be iterated until the disk thickness and effective temperatures both converge. 
However, we limit our procedure to one pass because experimentation has shown that the bulk 
of the effect is achieved after a single application. As discussed above, the inner 0.1 AU, where 
almost all of Fi n originates, is altered only slightly by the reprocessing calculation. 
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Fig. 2. — Shapes of the r = 0.03 surfaces of the inner, outbursting parts of the three disk models, 
showing the hot “volcanos” surrounding the central star (black dot). Scales on the horizontal and 
vertical axes are the same, so the disks are shown in their true proportions and the extent of 
reprocessing can be gauged. 
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Fig. 3. — Disk temperature at the r = 0.03 surface versus radius in the B1 model, with reprocessing 
(solid line) and without (dots). Also shown is the temperature at the base of the envelope for each 
radius (dashes). 




- 10 - 


2.4. Central Star 

Since the innermost radius r included in the disk models is 3 i? 0 , we must choose temperature 
and density distributions for r < 3i2 0 . If the central star is a T Tauri star of Solar mass and 
luminosity, its surface temperature and radius are between 3 000 and 4000 K and about 3ii 0 , 
respectively, implying that during outburst, the star’s own luminosity is negligible compared with 
that of the inner disk. However, the structure of the region where the disk meets the star is 
not known. Does the disk engulf the star with hotter material in outburst, as suggested by the 
two-dimensional radiation hydrodynamic calculations of Kley & Lin (1996)? The disk model does 
not account accurately for the interface between disk and star. It assumes Keplerian rotation even 
at r = 3 Rq, whereas at this radius the radial pressure gradient is large, suggesting some departure 
from Keplerian rotation. At r = 35 0 , the model disks’ half- thicknesses are about 1 i? 0 . 

As a test of the importance of the innermost 3i? 0 to the integrated spectrum during outburst, 
two spectra were calculated as described in § 3.2. For the first, the central 3i? 0 was occupied 
by a sphere of uniform surface temperature 3400 K. For the second, the vertical structure of the 
disk model at 3 i? 0 was copied to smaller radii. The models’ effective temperatures at 3 i? 0 are 
5-6000 K, while their midplane temperatures are 1.0-1. 2 x 10 5 K. The resulting spectra differ by 
at most 0.016 decade in flux over the entire wavelength range. In the remainder of this paper, the 
center of the disk is occupied by a star of temperature 3400 K and radius 3i? 0 . 


2.5. Envelope 

FU Orionis stars show flux excesses at infra-red wavelengths A ^ 10 fim when compared both 
with normal giant stars and with equilibrium dusty disks. In V1057 Cygni, the excess diminished 
after outburst in step with flux at 5-band, suggesting the 10 /im flux may be due to absorption 
and re-emission of radiation from near the central object by material out of the plane of the disk, 
such as a circumstellar envelope (Kenyon &; Hartmann 1991). As we have seen in § 2.3, the excess 
cannot be fully explained by disk reprocessing, because reprocessing does not heat the disk surface 
outside about 10 AU radius. For a viewer on the surface outside this radius, the outburst region 
is hidden by the curve of the disk. 

Our simplified model of the envelope is a layer of uniform thickness, touching the top of the 
disk, and with a central hole exposing the inner disk (Figure 1). This structure for the envelope is 
suggested by the results of two-dimensional collapse calculations, in which the density structure of 
matter infalling onto a disk is closer to plane-parallel than to spherical geometry, and most of the 
optical thickness is just above the surface of the disk (Yorke, Bodenheimer, & Laughlin 1993). Our 
envelope is made from material with the same opacity function as that in the disk (§ 3.1). The 
temperature Tj in the envelope a distance d from the center of the system is calculated assuming 
the envelope dust is in approximate radiative equilibrium with the object’s hot inner region. That 
is, the rate of energy loss by radiation is balanced by the rate at which radiation from the hot 
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region falls on the dust: 


47rcrT d 4 



(3) 


where Lu is the “illumination luminosity”, the portion of the luminosity of the outbursting part 
of the disk which is intercepted by the envelope. This temperature distribution is correct for an 
optically thin medium. For an optically thick medium, equation (3) is also approximately correct. 
This is indicated by the radiative transfer calculations of Hartmann, Kenyon & Calvet (1993) 
for spherical optically thick envelopes around Herbig Ae/Be stars, and of Kenyon, Calvet & 
Hartmann (1993) for spherically symmetric protostellar envelopes in which density is proportional 
to d -3 / 2 . The variations of temperature with distance resulting from these calculations range from 
T ~ d -0,75 to T ~ d -0,4 . In non-spherically-symmetric models Yorke et al. (1993) find T ~ d -1 / 2 
also, in both optically thick and optically thin regions. For a thin disk, T ~ d -1 / 2 reduces to 
T ~ r~ 1 / 2 , which yields a flat A F\ spectrum (Adams, Lada & Shu 1988) as observed in FU Orionis 
objects at A ^ 5 p m. 


The parameters which fix the geometry of the envelope are its thickness Az parallel to the 
disk axis, the radius of the central hole and the extinction Ay NV through the thickness. At 
the inner edge of the envelope, Az is larger than the disk thickness. Density in the envelope is 
chosen so as to yield the specified extinction at temperatures between 125 K and 800 K. Opacities 
in this temperature regime are described in § 3.1. The thickness of the disk beneath the envelope 
is determined including reprocessing, but because the disk thickness is changed at most 10% by 
reprocessing, the effect on the spectrum is negligible. 

When the optical depth through the envelope is less than unity, the spectrum of the envelope 
appears as emission or absorption lines superposed on the spectrum of the disk. These lines, 
due to water ice and silicates, appear at wavelengths 3-30 /im. Such lines do not appear in 
spectra of FU Orionis and V1057 Cygni (Cohen 1980) (except during dimming events such as 
that experienced by V1057 Cygni in 1995 - Wooden et al. 1995; Wooden 1996). The lines are 
absent from the models’ spectra provided that at each radius the envelope has an optical depth 
greater than unity at wavelengths where it and the disk beneath it contribute significantly to the 
integrated SED. When this condition holds, the spectrum is also insensitive to both the density in 
the envelope and its variation with radius. The condition is satisfied when the extinction through 
the envelope is uniformly above Ay NV = 100 magnitudes. It is satisfied almost everywhere when 
density varies as the inverse square root of the disk radius, with Ay NV = 100 magnitudes at disk 
radius r — 1 AU. Such a density distribution is obtained when a spherical p ~ d“ 3 / 2 distribution 
is flattened vertically into a disk. The region where the condition is not satisfied in this case is 
outside 25 AU, where the envelope is now optically thin at wavelengths ^ 30 pm. Compared with 
the case where the condition is satisfied throughout the envelope, this results in a spectrum which 
falls off more rapidly as A increases past 30 ^m. Throughout the remainder of this paper, density 
in the envelope is held uniform and equal to the value which results in Ay NV — 100 magnitudes. 
Resulting densities lie between 10" 14 and 10~ 12 g cm“ 3 , for the parameters listed in Table 3. 
Corresponding mass infall rates can be estimated using the results of the calculations of Yorke et 
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al. (1993). Typical velocities perpendicular to the disk surface for the infalling material in their 
simulations are about 3 km s _1 . When the infalling material has a density of 10 -14 g cm” 3 over 
the inner 25 AU of the disk, this yields a mass infall rate of 10 “ 5 M0 yr -1 , similar to the rates 
used in the time-dependent disk models (Table 1). 

If the model FU Orionis objects were fully self-consistent, the envelope temperature as set 
by La would be continuous with the temperature of the disk hidden beneath it. However, as 
Figure 3 indicates, while the temperature calculated from reprocessing agrees quite well with the 
dust equilibrium temperature assigned at the same radius throughout most of the disk, there is 
significant departure outside r — 10 AU, where the curve of the disk hides the hot inner regions 
from points on the disk surface. Here we maintain illumination of the envelope since the radiation 
transfer calculations of Yorke et al. (1993) indicate that in the infalling material T ~ r” 1 / 2 , 
independent of the optical depth along the line of sight to the central source. 


3. Calculation of emergent radiation 
3.1. Wavelength-dependent opacities 

New opacities for use in the radiative transfer calculation were obtained from Alexander (1995), 
calculated as described in Alexander & Ferguson (1994) for the S92 composition (Seaton et al. 1994) 
with X = 0.70 and Z — 0.02. All known sources of opacity which are significant at temperatures 
between about 700 K and 12 500 K are included. Alexander convolved monochromatic opacities 
with a Gaussian filter of half-width equal to 2% of the inverse wavelength, because of the very 
broad range of wavelengths to be covered. The result is a table of opacities k^(R,T) averaged 
around 100 wavelengths A. R is equal to p/{T/ 10 6 K) 3 . The 100 wavelengths are evenly spaced in 
log A from 100 nm to 100 pm. Spectral resolution is thus A/AA — 14, sufficient to reveal individual 
strong lines. For example, Ha at wavelength 656 nm is represented by one opacity value on the 
line, and adjacent continuum values. The table has entries every 0.5 decade in R and every 
0.2 decade in T, except for temperatures from 800 K to 3300 K, which bracket rapid changes in 
opacity due to dust evaporation and ionization. Here there are entries every 0.1 decade in T. A 
few points required for the radiative transfer calculations fall outside the available R and T. In 
these cases, we use k\ at the same T and the nearest available R. 

Below 700 K, there is only sparse data available on the optical properties of the materials 
which condense as dust grains, and the opacities are less reliable. In the absence of better 
frequency-dependent low-temperature opacities, Alexander’s 800 K opacities are used in the 
present calculation down to a temperature of 125 K. These opacities include grains made of 
silicates, iron, amorphous carbon, and silicon carbide. Graphite is not included in the Alexander 
opacities, presumably because its importance in protostellar environments is controversial (Pollack 
et al. 1994). Yorke (1979) also provides opacities in this temperature range, but for silicate grains 
only. Neither Yorke’s data, nor the extrapolated Alexander data has temperature or density 
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dependence. The Yorke values are slightly lower than the total opacities of Alexander at most 
wavelengths. For temperatures of 125 K and below, water ice becomes an important source of 
opacity. We use the sum of Yorke’s (1979) opacities due to water ice and silicates. 


3.2. Spectra 


Radiation emerging from a distribution of density p(r, z) and temperature T(r,z) is calculated 
by integration of the equation of radiation transfer (Mihalas 1978) 


dh 

d.T X 


= h~S x 


( 4 ) 


along a set of lines of sight through the distribution. The numerical procedure is a modification of 
that described by Yorke (1986). The source function used is a black body, S\ = B\(T). The optical 
depth t\ is defined by dr\ = K\pdx, with wavelength-dependent opacity K\{p,T). Equation 4 is 
solved separately at each wavelength. Wavelengths used are those for which opacities are provided 
(§ 3.1). For a given line of sight, optical depth t\ = 10 is located by integrating from the observer 
into the medium. The outward integration, which determines the emergent specific intensity, 
then begins at this point. The specific intensity I\ is set equal to the black-body intensity B\ 
corresponding to the local temperature, and equation 4 is integrated back towards the observer. 
At each integration step, the step length Ax is chosen so that its optical depth A t\ = n\p Ax 
is at most 0.1, and also so that it is shorter than one third of the local spatial resolution of the 
outburst model. If this step would leave or enter the material making up the disk or envelope, 
the step length is shortened so as to end the step within 10~ 5 Rq of the boundary. This ensures 
physically thin layers have the correct optical depth, and abrupt opacity and temperature changes 
near boundaries are well-represented in the integration. As needed for finding the source function 
and opacity along the line of sight, density and temperature are linearly interpolated on the 
non-rectangular grid made by the series of vertical models. Opacity at each wavelength is found 
by interpolating in log space on the grid of densities and temperatures described in § 3.1. 


The integrated spectrum of the model is calculated by summing emergent spectra over 
regular square grids of lines-of-sight. Because temperature variation in the disk occurs over small 
spatial scales near the star, and larger scales further away, a nested, concentric set of line-of-sight 
grids is used. The innermost grid is 0.2 AU across, and successive grids are five times wider. 
The outermost extends to 125 AU radius. When viewed from directly over the pole, the disk 
is circularly symmetric, and the grids are one-dimensional, each with 51 lines-of-sight along the 
radius. For inclination angles 0 < i < 90°, the view of the disk is symmetric only about the plane 
containing the polar axis and the line-of-sight. In this case, each grid is two-dimensional, and 
consists of 51 lines-of-sight along the radius perpendicular to the symmetry axis and 101 along 
the diameter parallel to the symmetry axis. The resolution on the innermost grid is 0.004 AU, or 
0.85 -Rq. With the disk models described here, increasing the resolution beyond this level produces 
negligible change in the total spectrum. 



-14- 


The apparent magnitude of the model in the B spectral band is calculated by integrating 
in wavelength the product of the spectrum F\ and a filter transmission curve, and applying 
the calibration between integrated flux and magnitude of Colina & Bohlin (1994). The filter 
transmission curve used is that tabulated by Johnson (1965). 


3.3. Images 

The same radiative transfer calculation which yields the integrated spectrum of the object 
also yields a spatially-resolved image. Solution of the radiative transfer equation along a line of 
sight yields the spectrum of radiation emerging from the disk at one point. A flux is found by 
integrating in wavelength the product of this spectrum with a filter transmission curve. Fluxes 
from all the lines of sight are assembled on a grid, making an image of the object. 


3.4. Distances and Extinctions 

To match the absolute level of a physical model against that of an observed spectral energy 
distribution, we need the distance to the object and the amount of obscuring material along our 
line of sight. Measurements of these quantities for the three objects are collected in Table 2. 
The estimate of the distance to FU Orionis is taken from Murdin & Penston (1977), that of 
the distance to V1057 Cygni from Straizys et ai (1989). The distance to V1515 Cygni was 
estimated by Racine (1968). Extinctions to the three objects were assembled by BLHK. In adding 
extinction to the calculated spectra, we scale the wavelength dependence of interstellar absorption 
(Mathis 1990) so as to obtain the desired extinction at 0.55 //m, the central wavelength of the 
V band. 

Also in Table 2 are the luminosities of the objects and the corresponding models. The 
quantities listed are apparent B magnitude; de-reddened B-band luminosity; and de-reddened 
total luminosity, integrated from 0.380 (i m to 100 fim. B-band luminosities were calculated from 
the apparent magnitudes assuming Ry = Ay/E(B — V) = 3.1. Total luminosities were calculated 
by integrating the spectral energy distributions of Kenyon Sz Hartmann (1991), assuming for 


Table 2: Distances, Extinctions, and Current-Epoch Luminosities of FU Orionis Objects 


Star Distance Ay B Mag Lb/Lb.q L tot /Lo 



pc 

mag 

Data 

Model 

Data 

Model 

Data 

Model 

FU Orionis 

400 ± 60 

2.0 

10.6 


160 

160 

250 

260 

VI 057 Cygni 

550 ± 100 

3.1 

13.4 

13.79 

90 

50 

250 

200 

V1515 Cygni 

1000 ± 200 

2.8 

13.7 

13.59 

160 

170 

230 

290 
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V1515 Cygni a falloff at long wavelength with the same shape els that of V1057 Cygni. The models' 
B-band and total luminosities were calculated using the methods of § 3.2 and the parameters 
of Table 3. The discrepancy between the luminosities of the Cl model and V1057 Cygni was 
resolved by placing the Cl model at 500 pc rather than at 550 pc (see § 4.1). This lies within the 
uncertainty in the distance measurement. The apparent B magnitude of V1515 Cygni was found 
by extrapolating the light curve through a dimming which began in 1980 (Kenyon et al 1991). 
The S-band luminosity calculated for V1515 Cygni from the magnitude is thus not directly 
comparable with the total luminosity, which is calculated from a SED which includes the effects 
of the dimming. 


4. Results 
4.1. Spectral Fits 

Figure 4 shows our best matches to recent observed spectral energy distributions of the three 
stars. These are the data collated by Kenyon & Hartmann (1991), but without de-reddening. 
SED points are placed at the nominal wavelengths of the photometric bands, which in most cases 
are effective wavelengths for equal input energy at all wavelengths. Filter transmission curves are 
taken from Johnson (1965) for bands UBVRI and N\ from Bessell & Brett (1988) for JHKLM ; 
and from references under Kenyon Sz Hartmann (1991) for other bands. The models’ parameters 
are given in Table 3. 

Figure 4 also shows that the B1 model spectrum coincides with the V1515 Cygni data at 
R-band, and is an excellent fit at other wavelengths, if it is reddened by Ay = 3.2 mag instead 
of the nominal Ay = 2.8 mag. Some of the V1515 Cygni data were obtained while that object 
had not fully recovered from a dimming it experienced beginning in 1980. The dimming, which is 
thought to be a dust condensation event, was accompanied by an increase in the object’s reddening 
(Kenyon et al 1991). 

The data in Figure 4 include measurements obtained over about ten years, so a little caution 
is needed in matching them with models of any particular epoch. For example, points at 12, 25, 
60, and 100 jj,m were obtained with the Infra-Red Astronomy Satellite (IRAS) in 1983-84, while 
near-infrared data were obtained by Kenyon & Hartmann in 1989-90. The spectra are calculated 
from models at epoch 1993. 


4.2. Comparison with Other Spectral Fits 


In this subsection we compare the fits generated in § 4.1 with those obtained in previous 
attempts to model the spectral energy distributions of FU Orionis stars using accretion disks. 
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The models of Adams et al (1987) are constant mass-flux accretion disks, with excess IR 
emission arising in dusty shells, remnant wedges of the objects’ spherically-symmetric parent 
molecular cloud cores lying >1 — 10 AU from the central object. A spectrum, calculated assuming 
blackbody emission from the disk and envelope and grain opacities in the envelope, matches the 
spectrum of FU Orionis adequately at visible wavelengths, but rises from 10 to 30 /xm, in conflict 
with the data presented in Kenyon & Hartmann (1991) and plotted in Figure 4. 

The models of Kenyon, Hartmann & Hewett (1988) are also constant mass-flux accretion 
disks, but without a source of excess IR emission. Spectra are calculated by assigning each disk 
annulus the fluxes of a supergiant star of the same effective temperature, and summing over the 
disk. Annuli with temperatures below 3 500 K are assigned blackbody spectra. The resulting fits 
to FU Orionis and V1057 Cygni are adequate at wavelengths of 0.35-1 /xm, but too bright by 
almost a magnitude at 2 /xm, and too faint by over a magnitude at 20 /zm. 

In the flared-disk models of Kenyon & Hartmann (1991), a constant mass-flux disk reprocesses 
light from a central star. The surface of the disk is bent into a shape H ~ r 9 / 8 , where H is the 
disk thickness and r is the cylindrical radius. This makes reprocessing more effective than in a 
flatter, standard constant mass-flux disk. Spectra are calculated as in Kenyon et al. (1988). The 
fluxes at 3-30 /xm are higher than those from non-reprocessing disks, but in the wavelength region 
longward of 5 jum, which we match using a transition from disk to envelope, the best-fit flared-disk 
models of FU Orionis and V1057 Cygni are too bright at short wavelengths and too faint at 
long wavelengths. Kenyon & Hartmann (1991) also present spectra calculated from models of 
V1057 Cygni which include a spherically-symmetric dusty envelope and a steady accretion disk. 
The radial optical depth of the model envelopes at wavelengths 5-50 /im is close to unity. The 
spectra match the observations well, except that they show broad features at 3-30 /xm which do 
not appear in higher-resolution observations of FU Orionis and V1057 Cygni (see § 2.5). 

The SED fits in Figure 4 are more detailed in several ways than the three others just discussed. 
They are calculated from outbursting disk models in which the mass flux varies with radius; the 
models include reprocessing of disk emission by other parts of the disk; and the shape of the 
disk surface is consistent with the underlying accretion model. An envelope is included, based 
on results of hydrodynamic collapse simulations; and the SEDs are calculated from the models 
by radiative transfer with fairly complete frequency-dependent opacities. The more detailed 


Table 3: Parameters used in calculating the simulated spectra shown in Figure 4. 


Star 

Model 

Distance/pc 

i 

Lu/Lq 

Rht AU 

Az/AU 

FU Orionis 

Al 

400 

30° 

ii 

0.8 

0.8 

V1057 Cygni 

Cl 

500 

0° 

30 

0.1 

0.1 

V1515 Cygni 

B1 

1000 

0° 

33 

0.2 

0.4 
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Fig. 4. — Fits (solid lines) to current-epoch un-dereddened observed fluxes from Kenyon & 
Hartmann (1991) (points), for the three objects. Vertical error bars on the data points show 
the quoted uncertainties in fluxes. Horizontal error bars extend to the wavelengths where filter 
transmission drops to 50% of its peak value. Where the wavelength range is short or the flux error 
is small, error bars lie inside the points. Interstellar reddening is included in the model spectra. 
The B1 model spectrum is also reddened an additional 0.4 mag (dotted line), as a possible way to 
allow for V1515 Cygni’s sudden 1980 dimming. Parameters of the models are listed in Table 3. 
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calculation produces SEDs which match the observations closely, especially at wavelengths from 
5 to 100 /xm where the spectrum is determined largely by the envelope. The difference between 
observed and model fluxes is no more than about 0.25 mag at most wavelengths, with a few points 
off by up to 0.75 mag. 


4.3. Parameter sensitivity 

In the spectral fits of Figure 4, two sets of parameters vary the shape of the spectrum in two 
different wavelength regions. At visible and near-IR wavelengths, the flux is determined by the 
distribution of hot material, which is found only in the outbursting inner disk. At wavelengths 
longer than about 5 /zm, the luminosity is dominated by large surface areas instead of high 
temperatures, and the flux is set by the envelope. At intermediate wavelengths, inner disk, 
reprocessing by the outer disk, and envelope can all be important, depending on the parameters. 
As discussed below, the parameter set needed to fit a spectrum is not necessarily unique. 

The parameters which affect the disk portion of the spectral energy distribution are the 
viewing angle, i; whether reprocessing is included; the time-averaged mass flux M and mass of 
the central star, M*, to both of which the disk luminosity is proportional; and the interstellar 
extinction, Ay. Parameters which influence the envelope part of the SED are the luminosity La 
of the outbursting disk as seen from the envelope; the envelope extinction Ay NV , thickness A z, 
and central hole radius Rh\ and the viewing angle i. Several of these parameters are restricted 
by data, while others are unimportant to the spectra. Among the disk spectrum parameters, 

M* was set in the BLHK models to 1 M 0 for all three objects. This is based on a pre-outburst 
spectrum of V1057 Cygni indicating it was a T Tauri star and hence roughly of Solar mass, and 
space density arguments which suggest FU Orionis objects cannot be precursors of high-mass stars 
(Herbig 1977). The mass fluxes were chosen by BLHK so the light curves of the model disks would 
match the observed time variation of B-band light (Table 1). Interstellar extinction is not treated 
as a parameter; we use the values discussed in § 3.4. The angle at which we view V1057 Cygni is 
restricted by rotational line broadening measurements to i < 30° if the central mass is close to 
IMq, while the inclination of FU Orionis must be 25° < i < 70° (Kenyon et al. 1988). We choose 
to view the models pole-on ( i — 0°), except for FU Orionis, for which we assume an inclination 
of 30°. 

Among the envelope spectrum parameters, the extinction Ay NV is set to 100 magnitudes to 
ensure the envelope is optically thick at all wavelengths. Since the envelope has a temperature 
different from that of the underlying disk, an optically thin envelope would yield strong emission 
or absorption lines due to water ice and silicates at wavelengths 3-30 /xm. These are not observed 
in FU Orionis nor in V1057 Cygni (Cohen 1980; Wooden 1996). Kenyon & Hartmann (1991) 
used an envelope extinction of 50 magnitudes. However, the inner parts of our envelope are 
too hot for ice grains to exist, and consequently have reduced opacity. For an envelope with 
Ay NV = 50 mag, where ice is solid the optical thickness is high, but wherever ice has sublimated 
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0 1 0 1 
Log 10 [ Wavelength X / /xm ] 


Fig. 5. — Sensitivity of the spectrum of the B1 model to three disk parameters (left) and three 
envelope parameters (right). The model is viewed (a) 0, 20, 40, and 60° from pole-on (top to 
bottom); (b) with (upper) and without (lower) disk reprocessing (the envelope has been removed); 
(c) with disk luminosity scaled up (top) and down (bottom) by factors of two; (d) with luminosity 
Ln illuminating the envelope scaled by factors of 3 (top) and 1/3 (bottom) from its best-fit value; 
(e) with the radius of the hole in the envelope varied likewise (a larger hole leads to a higher flux 
at 3 /im); (f) with the thickness of the envelope varied likewise (a thicker envelope leads to a lower 
flux at 3 /im). See § 4.3. 
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the optical thickness is less than unity at wavelengths A > 4 /im. Increasing Ay NV beyond our 
chosen value of 100 mag has little effect on the spectra. Another envelope spectrum parameter, 
the luminosity L t i illuminating the envelope, is expected to be less than the total luminosity of the 
outburst region in the disk. Collapse calculations suggest that reasonable values for the envelope 
thickness A z may lie in the range 0.1 — 10 AU (Yorke et ai 1993). The hole radius R h must be at 
least 5 Rq if the hole is to expose the peak surface temperature in the outbursting disk, and must 
be 4Oi?0 « 0.2 AU if the entire outbursting region is to be visible. 

Variation of the models’ spectra over the remaining parameter space is illustrated in Figure 5, 
which sets out spectra calculated for the B1 model 

(a) at viewing angles i = 0, 20, 40, and 60°: Increasing the viewing angle reduces the 
projected area of the disk, hence decreases the flux at all wavelengths. Beyond a critical 
angle, the inner edge of the envelope hides the central outbursting region, and the flux at 
visible wavelengths plummets. 

(b) with and without reprocessing: The envelope is removed so the spectra reflect disk 
temperatures only. Reprocessing yields a spectrum which is almost flat at wavelengths 
of 3-10 fi m, whereas a standard constant mass-flux accretion disk with T ~ r -3 / 4 has a 
long- wavelength spectrum A F\ ~ A~ 4//3 . 

(c) with values of MM* equal to one half, once, and twice the standard value: 

The integrated disk luminosity is proportional to MM*. To scale luminosity L, we scale 
temperatures throughout the disk sis L 1 / 4 (Lin &; Papaloizou 1985). Note that this scaling 
does not accurately represent changes in spectral lines. Since the envelope illuminating 
luminosity Ln is not adjusted, the spectrum is unaffected at wavelengths longer than a 
few /xm. 

(d) with standard, one third, and triple luminosities illuminating the envelope: The 

level of the flat spectrum is determined by the luminosity illuminating the envelope. A high 
Ln means a high starting temperature in the envelope’s T ~ r" 1 / 2 , and so a high flux at 
given wavelength. 

(e) with holes in the envelope of the standard radius, one third, and three times 
standard: Increasing the radius of the hole exposes more of the disk material heated by 
reprocessing to temperatures of 1000-1500 K. Since the envelope at the same radius is cooler, 
exposing more of the disk raises the flux at 2 3 /im. 

(f) with envelopes of standard, one third, and triple thicknesses A z: Increasing 
the thickness of the envelope places the envelope’s upper, visible surface further from the 
illuminating source. This reduces the flux from the envelope at wavelength 5 /im, and total 
flux at this wavelength falls to near that of the disk alone. In addition, unit optical depth in 
the thicker envelope spans a larger range of temperatures, resulting in stronger absorption 
lines which are not observed in these objects. 
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The mapping between models and spectra is not one-to-one. There exist combinations of the 
parameters for which the calculated spectra are almost identical. For example, changing the 
distance from which we view the model multiplies the spectrum by a factor independent of 
wavelength, shifting it vertically on our logarithmic plots. This is almost the same as the effect of 
changing the angle from which we view the model, provided the viewing angle is not so large that 
the hot central region of the disk is obscured. Changing the mass inflow rate M, or equivalently 
the mass of the central star, M», scales the luminosity of the disk and so likewise shifts the 
continuous spectrum vertically. However, reasonable ranges in any of these parameters require 
only small changes in the accretion disk models. As the Cl model set up by BLHK is slightly less 
luminous than V1057 Cygni, we choose to place the model slightly closer than the star (§ 4.1). 
We could instead have increased its luminosity using the M, M scaling described above under 
item (c). Since scaling luminosity scales the peak surface temperature of the disk, it moves the 
peak wavelength of the un-reddened spectrum. However, once several magnitudes of interstellar 
reddening are applied, the peak in the \F\ plot always lies near the H- band, effective wavelength 
1.63 /im, as fixed by the reddening curve. 

The effects on the spectrum of the envelope hole radius Rh and thickness A z are similar and 
small, and they occur in overlapping wavelength ranges. The hole radius would have almost no 
effect on the spectrum if our model envelope were fully consistent with the reprocessing applied 
to the model disk (§ 2.3), in which case the disk temperature and envelope temperature would be 
nearly continuous. Thus the precise assumed geometrical structure of the envelope is not critical. 
For good fits to the observed spectra, the important requirements on the envelope are (1) Ay NV 
must be large, and (2) the temperature distribution must be r -1 / 2 , as in equation (1). 

The luminosity illuminating the envelope, Lu, is not degenerate with other parameters; 
the observed SEDs thus constrain the value of Lu for each object. The best-fit values listed 
in Table 3 are much less than the objects’ total luminosities in Table 2. This occurs because 
equation (3) assumes isotropic radiation from the central region, whereas in fact the bulk of the 
total luminosity is emitted from deep inside the volcano, towards the pole, and plays no part in 
heating the envelope. Compared with the region inside the volcano, the portion of the inner disk 
which provides the chief illumination for the envelope has about twice the surface area and half 
the surface temperature (figures 2 and 3). Application of the Stefan-Boltzmann law thus suggests 
the envelope will be illuminated by a luminosity around one eighth that of the total for the object. 
According to Table 2 and Table 3, the envelope intercepts between one sixth and one twenty-fifth 
of the total energy emitted by the object. 

More evidence that the disk models are compatible with the required Lu comes from close 
inspection of the disk surface shapes in Figure 2. The ratio of Lu (Table 3) to the total B-band 
luminosity of the object (which arises almost entirely in the outbursting region; Table 2) is highest 
for V1057 Cygni. The corresponding model, Cl, has the smallest rim to its volcano, exposing hot 
material to the largest range of radii in the outer disk. The next-highest Lu/Lb is exhibited by 
V1515 Cygni (Bl). Compared with FU Orionis (Al), its rim is smaller, its outer disk surface lies 
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higher, and its peak surface temperature occurs higher above the midplane. The three models 1 
ranking in openness of their central volcanos is the same as the ranking of the objects in the 
relative luminosities Luj L b needed to explain the observed envelope spectra. 

The models were truncated at a radius of 125 AU because the observed fluxes decrease with 
wavelength at A ^ 50 /im, corresponding to an envelope temperature of about 60 K. The disk 
surface temperature at the same radius is 10 K (Figure 3), comparable to that expected in the 
surrounding molecular cloud. The accretion disk model then breaks down because viscous heating 
is no longer the main energy input. Furthermore, radiative collapse calculations by Yorke et 
ai (1995) suggest that the disk plus envelope becomes optically thin near this radius. Thus the 
spectra would fall off beyond this wavelength even if the models were not truncated. 


4.4. Simulated 5-band images 

Figure 6 shows the appearance of the inner 0.3 AU of the disk in the B1 model. This is 
the model whose spectrum and light curve match that of V1515 Cygni, but here its envelope is 
removed and it is viewed 30° from pole-on, through a 5-band filter. Lighter areas of the images 
are brighter; the density of grey scales with the logarithm of the flux. The series of images is 
calculated from a series of B1 models spaced every 25 years through an outburst cycle. The series 
covers quiescence while material accumulates and slowly migrates in through the inner disk (1243, 
1543, 1943), through onset, when the temperature at the inner edge of the disk rises high enough 
for hydrogen to be ionized, through the outward propagation of the ionization front (1968, 1993, 
2093), and as the outburst declines (2168) and a new cycle begins (2193). Note the bright inner 
slope of the volcano, an asymmetric white crescent at the centers of the 1968, 1993, and 2093 
images. Also, note the increased reprocessing in the outer disk during outburst. Reprocessing 
of illumination from a central star has been added to the disk models when they are quiescent. 

In outburst, the disk dominates total light. In quiescence, the disk is fainter than the star. The 
radius of the assumed star is 35©, and its temperature is 3400 K. The images were computed as 
described in § 3.3. Though spatial scales this small in protostellar accretion disks are unlikely to 
be directly imaged in the near future, occultation measurements have already been used to resolve 
spatial scales of about 1 AU in the disk around T Tauri (Simon et ai 1996). 


5. Discussion and Conclusions 

In this paper, we have calculated SEDs from the outbursting accretion disk models of BLHK, 
using frequency-dependent radiative transfer and including disk-disk reprocessing. Previous 
workers have calculated SEDs for constant mass-flux disks by summing blackbody emission or 
stellar spectra over a f = 2/3 surface. Our spectra also provide a good match for the observed 
SEDs of FU Orionis, V1057 Cygni, and V1515 Cygni. In our model, the flux at wavelengths 
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Fig. 6. — Time series showing the progress of an outburst cycle in the model V1515 Cygni (Bl) 
disk, as seen through a 5-band filter. All frames cover 0.3 AU in radius, and share one flux scale 
which spans 18 decades. The brightest regions represent a 5-band flux through the disk surface 
of 1 x 10 11 erg cm~ 2 s -1 . The disk is tipped away from the viewer by 30° from pole-on. The Bl 
model SED shown in Figure 4 was calculated using the 1993 disk model in the fifth frame. A full 
outburst cycle in this model lasts 1150 years. 
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0.35-2 /x m is provided by the disk inside 0.25 AU. Flux at 2-5 /xm is provided partly by the disk 
inside 0.25 AU, partly by reprocessing in the outer disk of luminosity from the inner disk, and 
partly by a flattened envelope in approximate radiative equilibrium with a central luminosity L t i. 
Flux at 5-100 ^m is provided by the envelope. 

Our model FU Orionis objects have some remaining weaknesses. Since the luminosities 
needed to match the envelope spectra differ somewhat from those available in the models for 
reprocessing by the envelope, it is likely the shape of the volcano is not accurately rendered by 
the disk models (§§ 2.3, 2.4). However, among the three objects, a higher required envelope 
illumination luminosity does correspond to a more open model volcano (§ 4.3). Another hint of 
mild inaccuracies in the disk models comes from Figure 4, where both Al and B1 show excesses 
at 0.55 /xm over the corresponding measured fluxes. While the B1 excess may be due to the 1980 
dimming event in V1515 Cygni (§ 4.1), both excesses might also be removed if the steep shapes of 
the volcanos shown in Figure 2 were relaxed, reducing the small degree of reprocessing between 
the volcanos’ opposite walls. Whether this relaxation is realistic could be demonstrated using 
two-dimensional hydrodynamic simulations of the outburst region, allowing radial derivatives of 
the physical quantities to act on the structure. The volcanos would be flatter also if the viscosity 
parameter a in the outburst region were larger than the value of 10“ 3 used in the BLHK models. 
The disk models further omit the back reaction of self-illumination on the propagation of the 
outburst, which is expected to be a second-order effect. Again in Figure 4, the spectrum of the Cl 
model is about 0.1 decade too faint at 3 /xm, and 0.1 decade too bright at 5 /xm. If the model 
were more luminous or its volcano were more open, so that reprocessing were stronger, this mis-fit 
might be corrected by increasing the radius of the hole in the envelope. 

The model envelope is another area where further work is needed. Our model is an 
over-simplification chosen for the small number of parameters it introduces. Model envelopes 
could be generated using a two-dimensional radiation hydrodynamic collapse calculation, but the 
exploration of parameter space needed to find envelopes appropriate for the spectra would make 
this a large undertaking. However, reconciliation of the illumination luminosities with the central 
luminosities of the objects does require more realistic vertical density gradients in the envelope, 
such as would be obtained from collapse calculations. 

The disk models assume a constant, low viscosity parameter a = I0 -4 outside the outburst 
region. As discussed by BLHK, this small viscosity results in a disk which is self-gravitating for 
the given input mass fluxes at radii beyond about 1 AU, yet we extend the disks to radii of 100 AU 
and more. This clear inconsistency is perhaps justified because outside a radius Rh (0.8 AU in 
the largest case, the Al model), the disk is invisible, hidden beneath the envelope. The thickness 
of the outer disk does affect the spectrum via the envelope temperature, which is a function of 
distance from the central object (equation 3). But this distance is only a weak function of the disk 
thickness, so any effect on the spectrum is small. The other requirement on the outer disk by the 
outburst model is that the outer disk supply mass at a near-constant rate; the mechanism of mass 
transport there might for example be gravitational instability, rather than that in the standard 
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Shakura & Sunyaev (1973) accretion disk picture. 

Although the process which triggers the rapid rise-time outbursts of FU Orionis and 
V1057 Cygni is still not understood, work in this paper lends solid additional support to the 
hypothesis that a thermal accretion event in a protostellar disk is responsible for the outburst. 
The spectral fitting also shows that an optically thick dusty envelope is an important component 
of the system, required to fit the 5-100 fim region. The presence of the envelope underscores the 
relative youth of these objects. The calculations demonstrate that reprocessing of disk light by the 
disk itself can produce a surface temperature distribution T ~ r -1 / 2 , which would result in a flat 
spectrum in the (A, A F\) diagram even without the envelope. However, this effect is cut off when 
the disk surface becomes concave towards the midplane. This occurs at a radius of about 10 AU, 
where the temperature at the surface f = 2/3 reaches the ice condensation point. As a result, the 
envelope is still needed to explain the flat spectrum over the full observed range of wavelengths. 

Finally, we note that many more details of the mechanism powering FU Orionis outbursts 
can be extracted from time-resolved spectroscopic observations. Though available data are sparse, 
it is already known that various spectral features appear and disappear during the evolution of 
the individual objects. Useful in the future would be spectra from 0.1-100 /im, spaced in time to 
resolve rapid changes in the light curve. For example, BBW76 is a recently-discovered FU Orionis 
object which is currently fading (Eisloffel, Hessman &; Mundt 1990; Reipurth 1991). Once a light 
curve is available covering a sufficient fraction of an outburst, a disk model can be constructed for 
this object. Comparison of SEDs calculated from the model against those observed over time will 
provide an even more rigorous test of the validity of the outbursting accretion disk model. 

These calculations were possible because David Alexander graciously calculated for us 
the frequency-dependent opacities. Scott Kenyon was quick to answer our questions. We also 
learned from discussions with Pat Cassen, Martin Cohen, Doug Lin and Diane Wooden, and the 
recommendations of an anonymous referee. The work was supported by grants NAGW-3408 and 
NAGW-4456 from the NASA Origins of Solar Systems Program, and by a special NASA theory 
program which supports a joint Center for Star Formation Studies at NASA-Ames Research 
Center, University of California, Berkeley, and University of California, Santa Cruz. The research 
made use of NASA’s Astrophysics Data System Abstract Service, and the Simbad database, 
operated at CDS, Strasbourg, France. 
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ABSTRACT 

A model for the magnetic dynamo process in accretion discs is proposed that does 
not involve underlying turbulent motions. The key ingredient of the model is the 
magnetic buoyancy. Combined with Coriolis twist, the buoyant motions result in a 
backcoupling between azimuthal and radial components of the field. It is shown that 
for reasonable values of the free parameters of the model all three components of the 
magnetic field reach a stable saturation state consistent with the basic assumptions 
underlying the model. 

Key words: accretion, accretion discs - MHD. 


1 INTRODUCTION 

To amplify and maintain an originally weak magnetic field, 
the classical kinematic dynamo model (see e.g. Parker 1979, 
Section 19) employs both a regular large-scale and a chaotic 
small-scale (turbulent) velocity field. The large-scale velocity 
field is just rotation, whereas the turbulent velocities are 
usually associated with thermal convection. As a by-product 
of the dynamo action, a fluctuating small-scale magnetic field 
is created, whose amplitude is assumed to be much smaller 
than the characteristic intensity of the large-scale (mean) 
field. Another key assumption of the model is that both 
large-scale and turbulent flows are not affected by magnetic 
fields. Pudritz ( 1980), who calculated non-linear effects for a 
kinematic dynamo in an accretion disc, found, however, that 
fluctuating fields much more energetic than the mean field 
develop, strongly reducing the intensity of turbulent motions. 
As a result, the mean field cannot grow above a certain 
critical value, at which the turbulence is shut down. The 
mean field was estimated to saturate at ft < 75 (where fi is the 
ratio PJP m of gas to mean field magnetic pressure), i.e. far 
away from the thermal equipartition value of ft = 1 suggested 
as a natural dynamo limit (Galeyev, Rosner & Vaiana 1979). 
Recently, more serious problems with the classical dynamo 
have been reported. A study of the Galactic dynamo by 
Kulsrud & Anderson (1992) revealed that rapidly growing 
fluctuating fields may substantially modify the turbulent flow 
before any significant amplification of the mean field is 
achieved. A similar pessimistic conclusion was reached in a 
general study of astrophysical kinematic dynamos by 
Vainshtein & Cattaneo (1992): fluctuating magnetic fields are 

* UCO/Lick Observatory Bulletin No. 1318. 


readily generated at tiny scales (where they reach equipar- 
tition values), while the mean field remains extremely weak. 

Objections against the kinematic dynamo have also been 
raised by Hawley & Balbus (1992), who note that the classi- 
cal model may not be applicable to accretion discs, as in that 
case even weak fields are able to significantly influence the 
fluid flow, thus invalidating the assumption on which the 
kinematic dynamo is based. This effect is due to the 
magnetorotational instability (hereafter MRI, a term coined 
by Goodman & Xu 1994), whose importance has been 
revealed in a series of papers by Balbus & Hawley (see 
Hawley, Gammie & Balbus 1995 and references therein). 
Yet another, long-known objection against the kinematic 
dynamo in accretion discs is that the origin of the turbulent 
velocity field remains unknown. The only commonly 
accepted source of turbulent motions in accretion discs, the 
convective instability, is usually effective in limited regions of 
a disc only (see e.g. Stella & Rosner 1984), and thus it cannot 
solve the problem. 

The idea that magnetic fields could be amplified without 
any assumed (or underlying) turbulence is not new. Its origin 
goes back to the early 1970s, when it was noted that the 
magnetic buoyancy, which is capable of generating turbulent 
motions by itself, might be an important factor in the dynamo 
process (Parker 1971). Meyer & Meyer-Hofmeister (1983) 
were the first to employ that idea for an estimate of the 
saturation state parameters in accretion discs. The effect of 
magnetic buoyancy on accretion discs was also discussed by 
Galeyev et al. (1979) and Stella & Rosner (1984). The first 
complete disc dynamo model entirely independent of any 
underlying turbulence was proposed by Tout & Pringle 
(1992). In that model, the Parker instability (Parker 1979, 
sections 13.4 and 13.5) due to the azimuthal field com- 
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poncnt, and the magnetorotational instability due to the 
vertical field component arc employed to close the dynamo 
loop. Almost simultaneously, a Galactic dynamo model 
based on magnetic buoyancy was proposed by Parker 
( 1992). Finally, in a recent stellar dynamo model of Tout & 
Pringle (1995) the dynamo feedback is provided by the 
differential stellar rotation combined with turnover flows’ 
generated by magnetic buoyancy. 

Based on the above discussion, we find that a dynamo able 
to amplify magnetic fields without any underlying chaotic 
small-scale motions is a very promising idea. In particular, 
under the conditions characteristic of an accretion disc 
interior (strong shear, strong vertical stratification due to 
variable gravitational acceleration in the direction perpendi- 
cular to the mid-plane) a dynamo powered by hydromagnetic 
instabilities that tap the energy directly from the orbital 
motion becomes a possible solution to all problems of the 
classical kinematic scenario. In the present paper we discuss 
a simple disc dynamo model based on processes outlined in a 
recent paper by Rdzyczka, Bodenheimer & Lin (1995), who 
performed a numerical study of buoyancy effects in accre- 
tion discs. Our approach is generally similar to that of Tout & 
Pringle (1992, 1995). For reasons explained in Section 2, 
however, we do not rely on the MR1 to close the dynamo 
loop, and the dynamo feedback is provided by magnetic 
buoyancy and helicity of the motions the buoyancy induces. 

We also differ from Tout & Pringle in the assumptions con- 
cerning the field loss from the disc, for which we assume a 
more general form. Our dynamo is described by a set of equ- 
ations which are introduced and discussed in Section 2. The 
solutions of the dynamo equations are presented in Section 
3. in Section 4 we summarize the results and discuss their 
implications for the structure and dynamics of an accretion 
disc. 

2 THE DYNAMO EQUATIONS 

To explain the physical basis of our model we shall begin 
with a brief summary of the simulations performed by 
Rdzyczka et al. ( 1 995), who have followed the evolution of a 
seed poloidal field in an accretion disc under the assumption 
of axial symmetry. They find that the azimuthal component 
of the field B $ grows due to shear until the thermal equipar- 
tition value of (5= 1 is approached. The growth of is 
stopped by formation of buoyant plumes with vigorous 
internal motions, leading eventually to partial ejection of the 
field from the disc, fragmentation of the seed and the 
establishment of a turbulent velocity field. Shortly afterwards 
a quasi-stationary state is reached, in which the poloidal field 
assumes a patchy structure, with the medium between 
patches being either more weakly magnetized than the 
patches themselves or entirely field-free. 

Rdzyczka et al. (1995) argue that the patches are likely to 
follow the same evolutionary pattern as the initial seed, and 
they point out that if the assumption of axial symmetry were 
relaxed, the Coriolis force acting on buoyant elements of disc 
gas would result in helical twisting of the azimuthal field lines 
and amplification of the poloidal field. An example of three- 
dimensional effects that could counteract the dynamo action 
instead of supporting it is given by I lawley et al. ( 1995), who 
find that weak azimuthal fields are prone to the MRI which 
results in magnetoturbulence. Therefore, it is conceivable 


that B 0 could saturate far below' the thermal equipartition 
value, before magnetic buoyancy becomes important. 
Hawley et al. (1995) point out, however, that as time 
advances the wavelength of the dominant azimuthal mode of 
the MRI becomes longer, and indeed their fig. 8 shows the 
field organizing itself into more coherent structures which 
may eventually evolve into nearly axisymmetric patches. 

The simulations of Rdzyczka et al. (1995) indicate that, 
because of the patchy structure of the poloidal field in the 
saturated state, the MRI would not significantly influence the 
vertical component of the field /?_, and that the generation of 
the radial component /? r from B : due to MRI would be rather 
inefficient. Thus we assume that the key link in the dynamo 
loop is not the B.~*B r transformation due to the MRI [as 
proposed by Tout & Pringle (1992)], but the B^ B r transfor- 
mation due to the helical twisting of B p lines. Further 
discussion of this issue is deferred to Section 4. Here we 
would only point out that the MRI is by no means unimpor- 
tant in our scenario, as it causes the magnetized region to 
spread radially, resulting in the expansion of an initially 
localized dynamo process both towards the centre of the disc 
and towards its outer edge. 

Both Rozyczka et al. (1995) and Hawley et al. ( 1 995) find 
that the field is likely to change its orientation on a length- 
scale which is small compared w ith the disc thickness. This 
finding corroborates another assumption made by Tout & 
Pringle (1992), according to which the main source of 
magnetic flux loss within the disc is reconnection. Since, 
however, according to Rdzyczka et al. (1995) and Hawley et 
al. (1995) all three field components may change sign on 
scales that arc small in at least one direction , we do not 
restrict reconnection to B z [as was done by Tout & Pringle 
( 1 992)), but allow for it also in the case of B ^ and B r 

Following Tout & Pringle (1992), we write our equations 
in terms of B 0y B. and B r , which should be regarded as local 
average values of the field components. As indicated in the 
preceding paragraph, we assume that every component of 
the field is destroyed by reconnection. The reconnection 
proceeds at a rate given by the ratio of the momentary value 
of the component and an appropriate time-scale r lct . (which 
may be different for different components). Additionally, and 
in accordance with Tout & Pringle (1992 ), B+ and B r are lost 
from the surface of the disc due to the undulating modes of 
the Parker instability (Parker 1 979, Section 1 3) at rates given 
by the ratios of their momentary values to a buoyant time- 
scale r,». 

We assume that the disc is nearly Keplerian. Thus, B 0 is 
created from B, due to shear at a rate 

(i) 

t it T 

where I is the local Keplerian orbital period. The same 
Parker instability that results in loss of B ^ and B r from the 
disc generates the vertical component of the field at a rate 


dr r ( » t p 

The gain terms given by equations ( 1 ) and (2) are the same as 
in Tout & Pringle ( 1 992). To complete the dynamo we have 
to define the rate at which the radial field is generated from 



B Q . Our key assumption is that the buoyant elements acquire 
an angular speed equal to the cpicyclic frequency k. As in a 
Keplerian disc k is equal to the local angular orbital velocity 
Q, we effectively assume that the azimuthal field lines within 
the element are rotated on a time-scale of the local orbital 
period, i.e. that B, is generated at a rate 


dJKJK 

d/ / 
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Admittedly, this may seem to be an overestimate. Our B, 
generation rate is not, however, greater than that found by 
Tout & Pringle ( 1992), and furthermore the saturation value 
of our dominant field component (B^) turns out to be only 
weakly dependent on this assumed rate (see below). 
Assembling all gain and loss terms for each field component 
we obtain the final set of equations 
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motions is given by 


II 
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(7) 


where II is the half-thickness of the disc given by the 
standard expression 


// = 


fe<\ 
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( 8 ) 


V% is the Alfven velocity calculated from B^ and rj is a para- 
meter in the range 2-5 (Tout & Pringle 1992). When we 
integrated the equations with the formula (7) for r P , it turned 
out that /^ approached 0.5 B^ in the saturation state, which 
prompted us to adopt a more general, heuristic expression 


r,.= f] 


II 


(9) 


where V\ is the Alfven velocity calculated from B r That 
modification, however, did not result in any significant 
changes in the saturation state. 

The characteristic time-scales for reconnection, r of 
the azimuthal and radial fields are defined as 


rt: 


dz 1_ 

r ~v\ F ' 


(TO) 


where dz is an average vertical size (thickness) of a 
magnetized patch, and T is a parameter related to the 
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magnetic Reynolds number, in the range of 0.02-0.1 (Tout & 
Pringle 1 992). In the following, an effective thickness Az will 
be consistently used instead of the ratio <$z/T. Following Pout 
& Pringle (1992) we assume that the vertical field can 
reconnect both in the <j> and r directions, on characteristic 
time-scales and rjj c , respectively. The effective reconnec- 
tion time-scale rf ec is then given by 

<^. c )~' = (t&V'+(tZ x Y HD 

Both and r*' t . are defined analogously to rfT in equation 
( 10), with Az replaced by characteristic length-scales A*£. and 
A;' c for the reconnection of the vertical field in the <j> and r 
directions, respectively. For an estimate of we modify the 
prescription of Tout & Pringle ( 1992) (their formulae 2.3.2 
and 2. 3. 4-2. 3.6) to obtain 


2 1 

3^2 



( 12 ) 


where Az replaces the MRI length-scale. Provided that the 
wavelength of the buoyant instability in the azimuthal direc- 
tion is shorter than the characteristic azimuthal extent of a 
magnetized patch, i.e. that the patch undulates along its 
azimuthal extent [which we expect to happen based on simu- 
lations reported by Rozyczka et al. (1995)], we may simply 
approximate A;* by Az. Finally, we relate Az to II by the 
formula 


Az = f //, 


(13) 


where e is a free parameter which, given the range of T in 
equation ( 10), could in principle vary in the range 0 < 6 ^ 1 0 
(at e= 10 the average geometrical thickness of a patch dz 
approaches //). However, independently of the field strength, 
there is always a small-scale structure clearly visible in 
meridional cross-sections of both the Hawley et al. (1995) 
and Rozyczka et al. (1995) models, indicating that dz in 
excess of ~0.1// is rather unlikely. On the other hand, very 
small e results in extremely large loss terms in equations 
( 1 4)— ( 16), practically inhibiting the dynamo effect. For these 
two reasons, we restrict e to vary between 0. 1 and 1 .0, i.e. we 
assume that the upper limit for the characteristic scale of 
magnetic structures generated in the disc is of the order of 

0.1//. " 

Taking V as a unit of time, and c\fy rip (where p is the 
local density) as a unit of magnetic field intensity, we can 
rewrite equations (4)-(6) in the form 
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where time (' and field components b f and b : are now r 
expressed in the units defined above. We decided not to 
incorporate the numerical factor of >/2 jt into the parameters 
e and 7] in order to preserve their straightforward relation to 
the disc half-thickness //. 
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3 SOLUTION OF THE DYNAMO 
EQUATIONS 

The set of ordinary differential equations consisting of 
equations ( 1 4 )— ( 1 6) was integrated numerically using the 
fifth-order Runge-Kutta method with adaptive stepsize 
control, as implemented in the Numerical Recipes routine 
odeint (Press et al. 1992), until an equilibrium (saturation) 
state was reached. To examine dependence of the equi- 
librium fields on initial conditions, numerical integrations 
were carried out for a wide variety of small initial fi elds. F ig. 
1 shows the results for initial fields, in units of cjAnp, of 
bj = b r = b. = 1 0 ~ 2 , 1 0 “ 6 and 1 0 " 1 % with r] = 2.5 and 6 = 0.1. 
The system moves to the same equilibrium in each case - 
that is, the saturation field does not depend on the strength of 
the initial field. Additionally, analytic equilibrium solutions 
to equations ( 14)-( 16) were found by setting time derivatives 
to zero and solving the resultant set of three algebraic non- 
linear equations. We found that the set has only one real, 
non-zero solution. For the (rj, e) pairs checked, the equi- 
librium values returned by the numerical integrator agree 
with the analytic solution to better than four significant digits. 

As a test of the stability of the equilibria, the integration 
was also run with initial conditions both larger and smaller 
than the equilibrium values by 1 per cent. The system came 
to the same equilibria as in the cases with small initial fields, 
i.e. the saturated state of the dynamo proved to be stable. 

To show the dependence of the saturation field on the two 
parameters i] and 6, numerical integrations were run to equi- 
librium over a grid of (tj, e) combinations. The grid is 1 1 
points on a side, evenly spaced in rj and in log, 0 e, and covers 
the region where 2^rj<5 and 0.1 <6 <1. The results are 
plotted in Fig. 2, and selected points are listed in Table 1. As 


can be seen in Fig. 2, the strongest component of the field at 
saturation is b ^ with b r smaller by a factor of about 2. For 
6=1, b : is only slightly smaller than b r . When e =0.1, b : is 
less than b r by 1-1.5 orders of magnitude. Also from inspec- 
tion of Fig. 2, dependence of the field components on r] is 
weak, because the largest loss terms are the reconnection 
ones, which are determined by 6. 

To demonstrate the importance of the reconnection terms, 
we removed them from the b ^ and b r equations (14) and (15), 
and repeated the integrations. The much smaller remaining 
loss terms due to the Parker instability led to large intensities 
of saturation fields, such that the magnetic pressure was 
comparable with or larger than the gas pressure over the 
(?/,€) region examined. For example, for rj = 2.5 and 6 = 0.3, 
the results were b $ - 1 .3, b r — 0.42 and b, = 0.27. 

Finally, we examined the sensitivity of the saturation fields 
to the Coriolis gain term in the b r equation ( 1 5) by increasing 
it by a factor of 2. As a result, the value of the dominant 
component b ^ increased by a factor of 1.25, while h r and b 
changed by factor of 1 .6- 1 .7. Effects of comparable magni- 
tude but opposite sign were observed in another integration, 
with the Coriolis term decreased by a factor of 2. 

4 SUMMARY AND DISCUSSION 

We demonstrated that our dynamo equations yield stable 
saturation states to which the amplification process 
converges after a few Keplerian periods. The saturated field 
strength is independent of the initial seed field strength. It 
does depend, however, on two parameters, one charac- 
terizing the growth time of undulations on azimuthal field 
lines t rj ). and the other the average thickness of magnetized 
patches in the meridional plane of the disc (f). .Slowly 
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Figure 1 . The convergence of the solutions of the dynamo equations to a unique saturation state. 
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Figure 2. The dependence of saturation field components on the parameters rj and r. Fields are in units of i\j4 np. 


Table 1 . Saturation field strengths versus 
Parker time-scale parameter t / and re- 
connection parameter f. Fields are in 
units of c\ J4np. Starting field com- 
ponents in each case were b 0 = b, = />. 
= 10 ,: 


V 

t 

h 

b r 

K 

2 

0.1 

0.091 

0.042 

0.0043 

2 

0.3 

0.23 

0.10 

0.029 

2 

1 

0.50 

0.20 

0.16 

2.5 

0.1 

0.093 

0.043 

0.0036 

2.5 

0.3 

0.24 

0.11 

0.025 

2.5 

1 

0.56 

0.23 

0.15 

5 

0.1 

0.097 

0.045 

0.0019 

5 

0.3 

0.27 

0.12 

0.015 

5 

1 

0.72 

0.31 

0.12 


growing undulations (large rj) lead to smaller surface loss 
rates of B 0 and B n and to larger saturation values of those 
two components. On the other hand, the rate of B : genera- 
tion is lower for larger yj , resulting in a smaller saturation 
value. The dependence of the saturation state on r] is rather 
weak, however; the effect of changing f is much more 
strongly pronounced. Smaller e enhances the reconnection 
loss and results in a general decrease of the saturated field 
strength. 

We showed that with reconnection loss terms removed 
from the B 0 and /^ equations (equations 14 and 15) the field 
saturates at high values corresponding to /?< 1. Let us note 
that at fi ~ 1 the Alfven velocity, w hich may be regarded as a 
reasonable approximation of buoyancy-generated velocities 
(Meyer & Meycr-Hofmeister 1983; Tout & Pringle 1992, 
1995) is comparable with the velocity of sound < s . The 
characteristic time-scale for the escape of the field from the 
disc is then Hji\ ~ 1 /Q. Moreover, for ft < 1 a rapid increase 
of dissipation due to shock waves (not accounted for in our 
model) is expected. Thus, /$= 1 seems to be a reasonable 


natural limit to field amplification. It may be noted that the 
same conclusion was reached in the case of a kinematic 
(convective) disc dynamo by Galeyev ct al. (1979), who 
pointed out that at /i < 1 the effects of magnetic tension 
would become strong enough to significantly suppress the 
convection. Our conclusion is that both the reconnection 
terms discussed in the present paragraph should be retained 
in the dynamo equations. 

In the saturation state the strongest component is B 0> with 
B r generally smaller by a factor of ~ 2, and B ; smaller by 
another factor of 3-10. Thus the structure of the field in the 
saturation state is entirely different from that obtained by 
Tout & Pringle ( 1 992), who found B ~ B 0 and B f * 0.1 5B r 
It is, however, consistent with the assumed structure of the 
field, which in a meridional cross-section of the disc may be 
described as an ensemble of radially elongated patches with 
/^constituting the main component of the poloidal field, and 
with the azimuthal field being still stronger. (Let us note that 
a patch of the poloidal field in the meridional plane corre- 
sponds to a bundle of toroidal field lines.) The basic 
uncertainty of our model concerns the Coriolis gain term in 
the /^ equation. We showed that an increased efficiency of B t 
generation results in a relative increase of B r with respect to 
B p, while a decreased efficiency results in a relative decrease 
of B : with respect to B r Since at the standard efficiency of 
Coriolis twist introduced in Section 2 B f seems to be already 
rather strong, and B rather weak compared with B 0 , we 
conclude that the Coriolis twist rate should not differ much 
from the assumed value of 1/7 if a saturated field with 
dominant and non-negligible B. is to be obtained. 

Its simplicity notwithstanding, our model produces 
entirely coherent results. First of all, reasonable saturated 
field strengths, somewhat below thermal equipartition, are 
obtained from reasonable assumptions concerning the range 
in which tj and e are allowed to vary. Secondly, the saturation 
values of all field components become consistently smaller 
with decreasing f, which means that fields w'ith significant 
small-scale structure are weak. This is exactly what one 
expects, as at a lower field intensity the magnetic tension is 
weak and the field lines may easily be wound up very tightly. 
Thirdly, the saturated field remains constant in time, as 
opposed to the 'Pout & Pringle (1992) model in which the 
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saturated field oscillates on an approximately orbital time- 
scale with an amplitude reaching - 60 per cent of the mean 
value. Fourthly, as mentioned by Rozyczka et al. ( 1 995), even 
though our model involves turbulent motions it is free from 
the basic problem of the turbulent kinematic dynamos that 
have been raised by Kulsrud & Anderson (1992) and 
Vainshtein & Cattaneo (1992). While in the kinematic 
dynamo the turbulence may prevent amplification of the 
mean field, in the buoyant dynamo it does not set in until the 
azimuthal field has become strong enough to induce buoyant 
motions. 

Tout <& Pringle (1992) obtain an oscillating field because 
of the shut-off of the MR1 at high /T, resulting in a vanishing 
gain term in the B r equation ( B r generation is resumed as 
soon as B : falls below a critical value due to continuous 
reconnection). For the shut-off to occur, however, the B, lines 
should extend throughout the whole thickness of the disc 
(connecting, possibly, to an external field). On the basis of 
numerical simulations reported by Hawley et al. (1995) and 
Rozyczka et al. ( 1995) we feel, however, that the internal disc 
field is more likely to be organized in rather thin and flat, 
radially extended patches. The effect of the MRI acting on a 
patch is just radial stretching, and there is no chance for a 
shut-off to occur. Thus we feel that the scenario proposed by 
Tout & Pringle (1992) is more appropriate for a disc 
permeated by an externally connected vertical field. If such a 
field is absent, i.e. if one is interested in amplification of 
purely internal disc fields, our scenario seems to be more 
appropriate. 

In conclusion we would like to stress that the simple 
model presented here is not so much meant to represent real 
astrophysical discs as to draw attention to hitherto 
unexplored possibilities. We merely demonstrated that a self- 
consistent scenario for field amplification and maintenance 
in accretion discs based on magnetic buoyancy and Coriolis 


effects is possible, and we believe that further research on 
this subject is highly desirable. 
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ABSTRACT 

We follow the evolution of a closed loop of relatively weak poloidal magnetic field, originally embed- 
ded somewhat above the midplane in a hydrostatic accretion disk. The equations of magnetohydro- 
dynamics are solved on a numerical grid in axisymmetric geometry. Viscous heating, radiative transfer, 
and the horizontal and vertical components of the gravity of the central star are taken into account. As 
the evolution proceeds, toroidal field is built up as a result of shear in the disk, and the field becomes 
buoyant as a result of interchange modes of the Parker instability. The effective wavelength of the 
buoyant instability, and its dependence on the strength of the initial field loop, are found to be consis- 
tent with a linear stability analysis. The buoyancy results in turbulent motions and expulsion of some 
field from the disk. Eventually, a saturation state is reached, in which the field assumes a patchy struc- 
ture, and the ratio of gas to magnetic pressure stabilizes in the range 1-5. Outward angular momentum 
transfer and an accompanying radial expansion of the magnetized region occur as a result of magnetic 
torques, and an equivalent a-viscosity parameter is estimated. The implications of these results on the 
generation of a magnetic dynamo in a disk are discussed. 

Subject headings: accretion, accretion disks — instabilities — MHD 


1. INTRODUCTION 

The observational evidence for magnetic fields in accre- 
tion disks is in most cases circumstantial. The only disk for 
which measurements of magnetic field strength have been 
made is the primitive solar nebula, with fields of 0. 1-1.0 G 
estimated from the residual magnetization of carbonaceous 
meteorites (Brecher 1972; Butler 1972). However, there are 
good reasons to believe that accretion disks are permeated 
by fields of either external or internal origin. In the first 
case, the field can be generated by the central object (Shu et 
al. 1994, and references therein) or dragged into the disk 
from the ambient medium (Lubow, Papaloizou, & Pringle 
1994, and references therein). In the second case, the disk 
itself generates and maintains the field, owing to the 
dynamo effect (Tout & Pringle 1992, and references therein). 
In this paper we concentrate on the latter possibility. 

While it was pointed out almost a quarter of a century 
ago that torques from internal fields could be responsible 
for the anomalously high effective viscosity of accretion 
disks (Shakura & Sunyaev 1973), the details of the actual 
processes are still not too well understood. An important 
advance was made in a series of papers by Balbus & Hawley 
(see Hawley, Gammie, & Balbus 1995, and references 
therein, hereafter the series is referred to as BGH). These 
papers deal with a magnetorotational (term coined by 
Goodman & Xu 1994) instability that is able to amplify 
originally weak disk fields of external origin by up to several 
orders of magnitude, and to induce ordered “ two-channel ” 
flow (Hawley & Balbus 1992) accompanied by an efficient 
angular momentum transfer. However, the more likely non- 
linear outcome of the instability was found to be a hydro- 
magnetic turbulence, for which the Shakura & Sunyaev 
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2 Warsaw University Observatory, Al. Ujazdowskie 4, PL-00478 
Warszawa, Poland; mnm hydra.astrouw.edu.pl 
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4 University of California Observatories/Lick Observatory, Board of 
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effective viscosity parameter a ranged from 0.02 to 0.7, 
depending on the initial configuration assumed for the field. 
The most important conclusion of BGH is that with the 
help of magnetic fields one can induce both turbulence and 
angular momentum transfer in hydrodynamically stable 
regions of accretion disks. There are, however, two impor- 
tant physical effects that are not included in the simulations 
presented by BGH: magnetic buoyancy and resistivity. As 
the authors write, “buoyancy effects will undoubtedly have 
an important effect on the ultimate nonlinear fate of the 
disk ” (Hawley & Balbus 1992). 

Hawley & Balbus (1992) also note that the classical kine- 
matic dynamo model (Parker 1979, § 19) may not be applic- 
able to the amplification of internal fields in accretion disks, 
since even weak fields are able to significantly influence the 
fluid flow, thus invalidating the assumption on which the 
kinematic dynamo is based. Another, long-known objection 
against the kinematic dynamo in accretion disks is that its 
operation requires an underlying turbulent velocity field 
maintained by an essentially unidentified factor. The only 
commonly accepted source of turbulent motions in accre- 
tion disks, the convective instability, is usually effective in 
only limited regions of a disk (see, e.g., Stella & Rosner 
1984), and thus it cannot solve the problem. The idea that 
the magnetic buoyancy might be an important factor in the 
dynamo process was formulated by Parker (1971). Meyer & 
Meyer-Hofmeister (1983) were the first to employ that idea 
for an estimate of the saturation state parameters in accre- 
tion disks. The effect of magnetic buoyancy on accretion 
disks was also discussed by Galeyev, Rosner, & Vaiana 
(1979) and Stella & Rosner (1984). The first complete disk 
dynamo model entirely independent of any assumed turbu- 
lence was proposed by Tout & Pringle (1992). In that 
model, the Parker instability (Parker 1979, §§ 13.4 and 13.5), 
arising from the azimuthal field component, and the magne- 
torotational instability, arising from the vertical field com- 
ponent, are employed to close the dynamo loop. Almost 
simultaneously, a Galactic dynamo model based on mag- 
netic buoyancy was proposed by Parker (1992). Finally, in a 
recent stellar dynamo model of Tout & Pringle (1995), the 
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dynamo feedback is provided by the differential stellar rota- 
tion combined with “turnover flows” generated by mag- 
netic buoyancy. 

Recently, the classical kinematic dynamo has also been 
questioned for reasons quite independent of the BGH con- 
clusions concerning significant dynamic effects of weak 
fields in accretion disks. A study of the Galactic dynamo by 
Kulsrud & Anderson (1992) revealed that rapidly growing 
fluctuating fields may substantially modify the turbulent 
flow before any significant amplification of the mean field is 
achieved. A similar pessimistic conclusion was reached in a 
general study of astrophysical kinematic dynamos by 
Vainshtein & Cattaneo (1992): fluctuating magnetic fields 
are readily generated at tiny scales (where they reach equi- 
partition values), while the mean field remains extremely 
weak. Again, under the conditions characteristic of an acc- 
retion disk interior (strong shear and strong vertical stratifi- 
cation owing to variable gravitational acceleration in the 
direction perpendicular to the midplane), a dynamo which 
is powered by hydromagnetic instabilities, and does not 
require any underlying turbulence to operate, becomes a 
possible solution to the problem. 

Based on the above discussion, we expect the magnetic 
buoyancy to play a very important role in accretion disks. 
The buoyancy can modify significantly the magneto- 
rotational instability investigated by BGH, and it can con- 
tribute in a significant (possibly even dominant) way to the 
processes that amplify and maintain the internal magnetic 
field of a disk. Last but not least, it may be a crucial factor 
responsible for the generation of the disk’s corona (Galeyev 
et al. 1979). The latter issue, however, will not be considered 
here. The aim of the present paper is to follow the develop- 
ment of the buoyant magnetic instability in a realistic accre- 
tion disk environment, taking into account as much physics 
as possible given limited numerical methods and computa- 
tional resources. In particular, we attempt to find answers 
to the following questions: 

1. What are the physical details of the processes that 
transform the primary orbital energy into the field and 
kinetic energies? 

2. Is the magnetic buoyancy indeed a vital factor limiting 
the growth of the field inside the disk, as suggested by 
Galeyev et al. (1979)? 

3. Can the magnetization initiated by a weak seed field of 
limited extent spread over large regions of the disk? 

4. How sensitive are the properties of the saturated state 
to assumed initial conditions (seed field strength and 
location)? 

5. Can the effective viscosity of a magnetic disk be esti- 
mated? 

To that end, we perform numerical simulations of the 
evolution of a seed field in an accretion disk under the 
assumptions of axial symmetry and ideal magnetohydro- 
dynamics. The technical details of the simulations are dis- 
cussed in § 2. The results are presented in § 3, which also 
contains a discussion of physical processes responsible for 
the evolution of the seed field and a linear stability analysis. 
The results are summarized and discussed in § 4, together 
with some prospects for future research. 

2. simulations: technical details 

Following the line of thought pioneered by Shakura & 
Sunyaev (1973), let us assume that the effective viscosity in 


accretion disks is of hydromagnetic origin (BGH; Tout & 
Pringle 1992 and references therein). Alternative viscosity 
sources, such as thermal convection (Lin & Papaloizou 
1980; Kley, Papaloizou, & Lin 1993) or gravitational insta- 
bility (Paczynski 1978; Larson 1984; Laughlin & Bodenhei- 
mer 1994) will not be considered. The ultimate aim of MHD 
simulations then would be to obtain a self-consistent model, 
in which radiative energy loss is balanced by heat gener- 
ation arising from orbital energy release induced by 
magnetic torques. At the present moment, neither 
computational resources nor physical understanding of the 
processes involved seem to be sufficient to tackle that 
problem; therefore, we feel that a gradual approach is more 
appropriate. Our strategy is to begin with simple models, 
allowing for selected physical processes and limited inter- 
actions between the disk and the magnetic field, and to 
extend them step by step toward greater complexity. Thus, 
the simulations presented below are not directly related to 
any class of real astrophysical disks. Rather, they should be 
regarded as laboratory experiments set up to investigate the 
response of magnetic fields to an environment which con- 
tains the most important ingredients of a real disk’s interior. 


2.1. Numerical Methods and General Assumptions 

The simulations have been performed with the help of the 
two-dimensional hydrocode originally described by 
Rozyczka (1985) and recently adapted to deal with prob- 
lems of ideal MHD, following the recipes of Stone Sc 
Norman (1992). The equations solved are the continuity 
equation (1), the energy equation (2), the momentum and 
angular momentum equations (3), and the induction equa- 
tions (4). The coordinates used are Eulerian cylindrical (z, r, 
</>), and rotational symmetry with respect to the axis r - 0 
as well as mirror symmetry with respect to the plane z = 0 
are assumed. For these coordinates and symmetries, the 
code successfully passed all tests suggested by Stone & 
Norman (1992). The basic variables are density p, internal 
energy density c, specific angular momentum J, the azi- 
muthal component of the magnetic field and velocity 
components v z and v n together with corresponding poloidal 
field components B z and B r in the z- and r- direction, respec- 
tively. In equations (l)-(4) B z , B r , and B 0 are the com- 
ponents of the vector B s while v z , v r , and J/r are represented 
by the vector v: 

dp 

^ + V(p®) = 0 , (1) 


- + v-m= -P,v«,-VF rad -0£, 



+ (» • V)r 


— VP. — p\Q> + — (V x B) 
4 n 

x B - e zr4> div Q , 


SB 

— = V x (i> x B ) , 


( 2 ) 

(3) 

(4) 


where P g in equations (2) and (3) is the gas^ pressure, related 
to p and e by the equation of state, Q and E are the 
viscous stress tensor and the symmetrized tensor of velocity 
derivatives, respectively, ^rad is the radiative energy flux, 
and the factor e zr4> , which is assumed to be equal to unity in 
equations for v z and v r , and to zero in the angular momen- 
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turn equation (see the discussion below), has been intro- 
duced to simplify the notation. 

According to the research strategy outlined above, we 
focus on the dynamical aspect of the disk-field and field- 
disk interactions, and we adopt a simplified scenario for 
heat generation in the disk. We confine ourselves to the 
framework of ideal MHD and assume that the resistivity of 
the medium in which the field evolves is negligibly small, i.e., 
that the physical dissipation of the field is negligible. Conse- 
quently, in our approach the heat generated by reconnec- 
tion and/or dissipation of the field is omitted from equation 
(2), and the field can only affect the thermal budget of the 
disk through the P g dV work arising from motions it excites. 
Instead of the magnetic dissipation, we introduce an 
assumed a-type viscous dissipation, and we include heat 
gain from this effect. To satisfy these assumptions, equation 
(2) includes terms allowing for heat gain or loss through 
radiation transfer and heat gain through viscous dissi- 
pation. We assume, however, that the angular momentum 
redistribution arises only from magnetic torques, so that the 
viscous term is omitted from the angular momentum equa- 
tion. This issue is further discussed in § 2.2. 

The hydrodynamical part of the code is based on a 
second-order method with monotonicity constraints on 
advected quantities. The ME1D part employs a descendant 
of the constrained transport method (Evans & Hawley 
1988) to evolve the poloidal components of the magnetic 
field, with the original algorithm supplemented by upwind- 
ing in Alfven wave characteristics (Stone & Norman 1992). 
With two exceptions the upwinding is applied to all terms in 
the MHD equations that describe the evolution of, or are a 
consequence of, magnetic fields. The exceptions are the 
magnetic pressure term in the momentum equations, which 
is treated in the same way as the gas pressure term, and the 
advective term in the azimuthal field equation, which is 
analogous to the advection terms for all scalar hydrody- 
namical quantities. The MHD scheme conserves the zero 
divergence of the magnetic field very well, but it suffers from 
numerical reconnection of the field on the scale of the grid 
spacing. While usually small in the regions of smooth flow, 
the latter effect becomes significant in a turbulent medium 
(see Hawley & Balbus 1992, and § 3.4 of the present paper). 
The energy of the field that has been numerically dissipated 
and/or reconnected is lost from the grid. The geometrical 
constraints disable cyclonic motions in the disk, resulting in 
a lack of feedback between the azimuthal and poloidal com- 
ponents of the field. Nevertheless, we are forced to adopt 
them for the simple reason that they make the problem 
tractable with the present-day numerical methods and 
resources. 

2.2. Input Physics 

The environment in which the fields evolve is an annulus 
of a steady state x-accretion disk orbiting a point gravity 
source and embedded in a uniform ambient medium. The 
central gravity source is assumed to be massive compared 
with the disk, and the self-gravity of the disk is neglected. 
Following Rozyczka, Bodenheimer, & Bell (1994), the 
ambient medium and the gas within the disk are treated as 
separate fluids. The interface between them is defined as an 
equidensity surface, corresponding to the minimum density 
Pi encountered in the initial disk model. At every time step, 
each grid cell is labeled as belonging to the ambient medium 
if its density is lower than 1.01 p f , and as belonging to the 


disk otherwise. No flow is allowed across the interface from 
the ambient medium into the disk, but it follows from the 
labeling procedure that a residual mass transfer from 
the disk into the ambient medium may occur. Effectively, 
the ambient medium serves only as a means of imposing a 
constant pressure boundary condition on the surface of the 
disk, whereby the surface itself is allowed to move almost 
freely across the grid. As the density of the ambient medium 
is reset every time step to p h the residual mass transfer from 
the disk into the ambient medium results in a global mass 
nonconservation. However, this effect is very small, and in 
no case was more than 1 % of the mass initially contained in 
the grid lost by the end of a simulation. The temperature of 
the ambient medium (T amb ) is chosen so that the ambient 
pressure roughly matches the lowest pressure encountered 
in the initial disk model. Both in the disk and in the ambient 
medium, the equation of state is taken to be that of an ideal 
gas. The ambient medium rotates around the central mass 
with rotational velocity equal to the local Keplerian veloc- 
ity everywhere, and no other motions are allowed in it. To 
meet the constant pressure condition, the ambient medium 
is kept free of any magnetic fields by a procedure discussed 
in § 3.4. 

The initial disk model is in thermal equilibrium, with heat 
sources from viscous dissipation balanced by radiative loss. 
The viscosity coefficient is obtained from the local values of 
p, P ff , and Q according to the standard a-prescription 


where 

cj = 7 ^ , (6) 

P 

and D = J/r 2 is the local angular velocity (see, for example, 
Ruden & Pollack 1991). Only the shear viscosity is taken 
into account, while the bulk viscosity coefficient is set to 
zero. Artificial viscosity is not applied anywhere in the disk, 
and the a-viscosity does not operate in the ambient medium 
nor in the thin layer (1-2 grid cells) immediately below the 
surface of the disk. To obtain the viscous heat generation 
rate that appears in the energy equation, all components of 
the viscosity tensor (Tassoul 1978) are calculated. The vis- 
cosity coefficient (eq. [5]) is used in the energy equation as 
well as in the r- and z-components of the momentum equa- 
tion, where it is found to have a stabilizing influence on the 
flow. However, as we already mentioned in § 2.1, viscous 
torques are consistently set to zero in the angular momen- 
tum equation. We do this to avoid viscous flow's of the type 
discovered by Rozyczka et al. (1994) that could interfere 
with (or even quench) the flows excited by Lorentz forces. 
At the inner and outer radial boundary of the grid (r = r nVin 
and r = r max , respectively) all variables are simply reset at 
every time step back to their initial values. This procedure 
in fact has very little influence on the results of the calcu- 
lations, because ail simulations are stopped before the mag- 
netized region of the grid approaches the boundaries. 

For the radiation transfer, Rosseland mean opacities 
based on solar composition are included in the manner 
described by Bodenheimer et al. (1990). The main sources at 
low temperatures (less than 3000 K) are molecules 
(Alexander 1975) and grains (Pollack, McKay, & Christof- 
ferson 1985). For the calculations reported in this paper, the 
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minimum optical depth across the grid in the 2 -direction is 
~ 10 2 . Most of this optical depth is associated with the disk, 
which is embedded in a semitransparent ambient medium 
( T amb ^ U Thus, one may assume that the gas and radiation 
temperatures are equal everywhere, and that the radiative 
transfer equation may be approximated by the flux-limited 
radiative diffusion equation. The latter is solved by the 
alternating-direction implicit method described by Boden- 
heimer et al. (1990). The radiation transfer solver requires 
that the temperature at the boundary of the grid, T b , be 
specified. In the ambient medium we set T h ~ T amb , whereas 
in the disk T h is set equal to the local disk temperature. 

2.3. Initial Disk Model 

For the simulations presented here we employ an 
ot-model of a protoplanetary disk around a solar-like star. 
However, our calculations are not meant to apply specifi- 
cally to a protoplanetary disk, and there is no reason for 
selecting that particular model, except for the availability of 
detailed structural data. As we indicated at the beginning of 
§ 2, our intention is to look for conclusions concerning the 
response of magnetic fields to an environment containing 
the most important ingredients of astrophysical disks in 
general (strong shear, variable gravitational acceleration 
and radiative energy transfer). 

For the construction of the initial model, the method 
developed by Rozyczka et al. (1994) is used. First, accurate 
one-dimensional vertical structure integrations for different 
radii between r min and r max are carried out for a = 0.05, 
M = 10 6 M q yr“ *, and a central mass of 1 M 0 (see Bell & 
Lin 1994 for details). The whole set of one-dimensional 
models is then transferred onto the two-dimensional grid by 
means of a linear interpolation and allowed to relax toward 
a stationary state without magnetic fields. During the relax- 
ation, viscous heating and two-dimensional radiative trans- 
fer are included, but the viscous torques are omitted from 
the angular momentum equation. The essence of the relax- 
ation is the damping process, in which v r and v z are multi- 
plied every time step by a factor which is initially set to 0.95 
and is subsequently increased up to unity. The relaxation is 
stopped when the total kinetic energy of z- and r-motions 
has stabilized. The relaxed models exhibit weak convective 
patterns with typical velocities of 10“ 2 km s *, much below 
the velocity of sound (0.3-3.0 km s " l ). 

Two initial disk models are obtained for annuli (r min , r max ) 
of (0.2, 0.5) and (1.0, 1.5) AU, hereafter referred to as the 
inner disk and the outer disk , respectively. The inner disk is 
simulated on a grid consisting of 60 x 200 points uniformly 
distributed in (z, r) and providing a resolution of % 2 x 10 10 
cm. For the outer disk, grids of 65 x 120 and 125 x 230 
points are used, providing resolutions of ^6 x 10 10 and 
^3 x 10 10 cm, respectively. 

2.4. Initial Field Generation 

In a relaxed model an initial loop of poloidal magnetic 
field is generated, hereafter referred to as the seed loop. The 
seed field is not generated by any realistic physical process, 
and both its configuration and its intensity are arbitrary. 
The sensitivity of simulations to the initial setup of the mag- 
netic field is discussed in § 3.2 and 3.3. The loop-generating 
procedure assumes a distribution of azimuthal currents j^,, 
for which the azimuthal component of the vector potential 

A * x) - f v j7 (71 


is found, where x and x are vectors in the (z, r)-plane, and it 
is assumed that the region with nonzero current density has 
a small radial extent compared to its distance from the disk 
center. The poloidal B-field components are then obtained 
from the general formula 

B = V x A , (8) 

which, in an axially symmetric case, reduces to 


}_ djrAf) 

r dr ’ 


(9) 


and 


B= — 



( 10 ) 


The numerically computed divergence of the poloidal field 
thus obtained is zero almost to machine accuracy, and it 
does not change noticeably during the evolution. 

The seed loop is subsequently normalized to the desired 
field strength B 0 , which we identify with the maximum field 
strength to be encountered in the loop. As all seed loops are 
flat (with a typical width-to-length ratio Az/Ar ~ | ; see Fig. 
\a), B 0 is a good measure of the initial radial component B r . 
We carry out simulations for seed loops with B 0 of 0.3 and 
3.0 G, hereafter referred to as the weak seed and the strong 
seed , respectively. The weak seed loop is generated in both 
the inner and the outer disk, and the initial value of [I (the 
ratio of thermal to magnetic pressure) within it approaches 
^3 x 10 4 and ~ 300, respectively. The strong seed is gener- 
ated in the inner disk only, corresponding to /? ~ 300. For 
reasons explained in § 3.1, all loops are generated above the 
midplane of the disk. Let us also remark that one may 



Fig. 1. — Strong seed in the inner disk shown at (a) t = 0.2 yr and ( h ) 
t — 0.4 yr. (a) Contours of constant density iri the (r, z)-plane (solid lines ) are 
separated by A log p — 0.34 with log p min = — 1 1.0. Contours of constant 
P , the ratio of gas pressure to magnetic pressure ( broken lines), are spaced 
logarithmically between 100 at the edge of the loop and ^ 5 inside the loop. 
Arrows : Velocity vectors with length proportional to speed and with 
Khm ~ 0.15 km s '. The horizontal feature cutting the loop into two 
nearly equal parts is the interface between regions with oppositely oriented 
azimuthal field. ( b ) Contours of constant specific angular momentum (solid 
lines) are separated by A log J = 0.01. Broken line: Contour of constant 
p = 3000. Arrows: Velocity vectors with F max = 0.26 km s" *. The distance 
unit on the horizontal and vertical scales is 2.76 x 10 u cm. The left and 
right dosed J-contours correspond to a local minimum and a local 
maximum of J , respectively. 
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expect the weak seed in the outer disk and the strong seed in 
the inner disk to evolve similarly, because they have the 
same /Lvalue. Our results reported in § 3.2 demonstrate 
that this is indeed the case. 

The insertion of the seed field into the disk destroys the 
stationary equilibrium achieved as a result of the relaxation 
procedure. To approximately restore it, the internal energy 
density is adjusted locally according to 


where P m is the magnetic pressure and e 0 is the original 
internal energy density. The relaxation procedure described 
in § 2.3 is then repeated, with the contribution from P m 
taken into account. During the relaxation, the B-field is not 
allowed to evolve, i.e., all time derivatives of all B- 
components are set equal to zero. In general, the internal 
energy adjustment brought about by equation (11) brings 
the model so close to a stationary state that little relaxation 
is needed. It must be stressed here that the relaxed model 
with seed is not in equilibrium in the sense that it begins to 
generate the B ^ component as soon as the field is allowed to 
evolve. The aim of the somewhat tedious relaxation pro- 
cedure is simply to exclude any spurious effects from artifi- 
cially excited motions in the initial disk model. 



time [yr] 

Fig. 2. — Growth of the azimuthal field in the buildup evolutionary 
phase of the weak seed in the inner disk. Diamonds: Maximum B ^ value 
found in the grid at a given time. Curve: value calculated for that time 

from the linear formula (eq. [ 1 2]). 


3. simulations: results 

As stated in § 2.1, throughout the evolution the field is 
assumed to be ideally frozen into the gas, i.e., the disk gas is 
assumed to be perfectly conducting. Since the temperature 
of our initial disk models is rather low (200 K < T < 2000 
K), this assumption may hold only marginally, especially in 
the outer disk (see Fig. 1 a of Stepinski, Reyes-Ruiz, & 
Vanhala 1993). The neglect of the finite resistivity of the disk 
medium comes partly from our research strategy and partly 
from the conviction that, at the modest resolution provided 
by our grids, the effects of physical and numerical diffusion 
would be difficult to separate. However, the finite resistivity 
of the disk medium should definitely be taken into account 
in future experiments. 


3.1. Early Evolution 

In all cases, the early evolution of a seed loop inserted 
into an accretion disk follows the same general track. In the 
first phase (the buildup phase), the B $ component is gener- 
ated by shear at a rate nearly constant in time and approx- 
imately given by 


dB ± 
dt 



( 12 ) 


The growth of B ^ in the buildup phase of the weak seed in 
the inner disk is illustrated in Figure 2, where the maximum 
azimuthal field B% found in the grid is plotted against time, 
together with predictions based on equation (12). For a 
given time, the predicted value of the azimuthal field is 
obtained from the momentary maximum value of the radial 
component B™ found in the grid and from ft corresponding 
to the momentary location of B The two sets of data are in 
good agreement, demonstrating that the resolution we are 
able to achieve within the seed loop is sufficient to properly 
describe the linear evolution of the field. For z > 1 yr, the 
effects of motions within the seed and the effects of the 


overall expansion of the seed become visible, and the 
maximum B ^ in the grid falls below the value calculated 
from equation (12). At the end of the buildup phase, loss of 
field owing to numerical reconnection also becomes notice- 
able (see discussion in § 3.4). 

Let us note in passing that, according to equation (12), 
the sign of B ^ is determined by the sign of B r , which means 
that the azimuthal field in the upper half of the loop is 
directed oppositely to the azimuthal field in the lower part 
of the loop. The fact that B $ reverses its polarity inside the 
loop motivated us to insert the seed loops away from the 
midplane of the disk, so that the oppositely oriented field 
patches would form in the grid, and the interaction between 
them could be taken into account in our simulations. 

As soon as the azimuthal field appears in the disk, a 
magnetic torque, corresponding to the azimuthal com- 
ponent of the Lorentz force 


1 1 d 

— B r - — rB 4 
An r dr 


(13) 


begins to act on the disk gas. Let us stress again that the last 
two equations are approximate and they become invalid 
when B z and/or vertical gradients of B $ and ft cannot be 
neglected. The Lorentz force L $ is negative at the inner edge 
of the loop and positive at the outer edge, where inner and 
outer identify the location of the loop's edge with respect to 
the disk center. As a result, the torque associated with L 0 
causes the angular momentum to drift from the inner into 
the outer edge of the loop. In response, the loop stretches 
radially (both inward and outward) while essentially main- 
taining its vertical extent (Fig. 1). It is worth stressing that 
angular momentum redistribution in the loop results in a 
net flow of mass toward the center of the disk, associated 
with an orbital energy release, and that the orbital energy of 
the disk is mainly transformed into the energy of the 
azimuthal field. The radial expansion of the loop is a mani- 
festation of the magnetorotational instability discussed by 
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BGH. However, because we use different initial and bound- 
ary conditions, our results cannot be directly compared to 
theirs. The one exception is the circular poloidal loop 
described by Hawley & Balbus (1991) in their § 3.5 and 
illustrated in their Figure 9. Indeed, in that case the loop 
evolves mainly by radial stretching, just as our loops do 
throughout the buildup phase. 

According to equation (13), the evolution of the disk 
should not depend on the orientation of the field in the seed 
loop, i.e., two seed loops differing only in the orientation of 
the initial field should undergo the same evolutionary 
changes. To check the consistency of the MHD code, the 
calculation illustrated in Figure 1 was repeated with the 
direction of the initial field reversed. No differences were 
found. 

3.2. Buoyant Phase 

Eventually, because of continuous B ^ growth, the value 
of P inside the magnetized region approaches unity, and at 
the upper surface of the loop a clear wavy pattern appears. 
At this moment the loop enters the buoyant evolutionary 
phase, in which plumes of magnetized gas begin to rise 
toward the surface of the disk. We tested the credibility of 
our simulations in that phase by following the evolution of 
a weak seed in the outer disk at resolutions of 6 x 10 10 and 
3 x 10 10 cm. We find that throughout both the buildup 
phase and the early stages of the buoyant phase the results 
are hardly distinguishable. The high- and standard- 
resolution models begin to differ only at more advanced 
stages of the buoyant phase (Fig. 3), when the plumes arrive 
at the surface of the disk. However, even then the differences 
are insignificant as far as the global parameters of the model 
(flow pattern, total magnetic energy, and total kinetic 
energy of meridional motions) are concerned. The basic 
conclusion from this experiment is that our standard grids 
are able to resolve the effective radial wavelength k b of the 
buoyant instability. 

In Figure 4 different seeds evolving in the inner disk are 
compared at the moment when the first buoyant plume 
reaches the surface of the disk. Figure 4a shows the 
standard strong seed (cf. Fig. 1), originally extending 
over Ar = 1.54 x 10 12 cm, whose center is located at 
Zq = 2.8 x 10 11 cm above the midplane of the disk and 
r 0 = 5.46 x 10 12 cm from the rotational axis. The standard 
weak seed with the same original dimensions and location is 
shown in Figure 4b. Figure 4c shows a weak seed with 
standard z 0 and r 0 but extending over Ar = 2.5 x 10 12 cm, 
whereas a weak seed with the standard r 0 and Ar but orig- 
inally located at z 0 = 3.47 x 10 1 1 cm is illustrated in Figure 
4 d. The simulations demonstrate that k b is a function of 
both the seed strength and the location of the unstable 
magnetized area in the disk. It can be seen that at a given r 0 , 
k b grows larger with increasing seed strength B 0 (compare 
Figs. 4 a and 4b), increasing distance r from the disk center 
(inspect Fig. 4b or 4c) or decreasing distance z from the 
midplane of the disk (compare Figs. 4b and 4 d). The first 
effect is very clear, while the remaining two are only detect- 
able after a closer inspection of Figure 4. It may also be 
noted that k b in the outer disk is much longer than any of 
the wavelengths excited in the inner disk. 

The excitation mechanism of the buoyant instability will 
be discussed in § 3.5. Here we would like to point out 
another important property of magnetic buoyancy, namely, 
that radiative energy transfer has a significant effect on its 





Fig. 3. Weak seed in the outer disk shown at the end of the buoyant 
phase (/ = 4.55 yr). (a) Standard resolution: Az = Ar - 6 x lO 10 cm. (b) 
High resolution: Az - Ar = 3 x 10 10 cm. Contours of constant density in 
the (r, z)-plane {almost horizontal, widely spaced lines) are separated by A 
log p = 0.25 with log /J min - - 1 1.25. Contours of constant magnetic pres- 
sure ( closely spaced and irregular lines ) are separated by A log P m = 0.3 
with (log P m ) min ~ -3.9. Arrows: Velocity vectors with V ma% - 0.5 km s' 1 . 
The distance units on the horizontal and vertical axes are 1.84 x 10 12 and 
10 1 2 cm, respectively. 

efficiency. It has been recognized for a long time that the 
magnetic buoyancy is aided by radiative heat flow into 
the rising elements. To explain this effect let us note that 
in the absence of radiative energy transfer the elements 
would cool adiabatically, thus weakening the effective 
Archimedes force (see Parker 1975; Stella & Rosner 1984). 
Here we are able to demonstrate this effect by comparing 
two simulations of the strong seed in the inner disk (Fig. 5). 
With all initial and boundary conditions identical, one of 
the simulations allows for the radiative transfer (Fig. 5a), 
while in the other the radiative transfer and the viscous heat 
generation are both shut off (Fig. 5b). In the nonradiative 
case, the buoyant plumes indeed grow much more slowly 
(while the wavelength of the buoyant instability does not 
seem to be affected to any significant degree). Unfor- 
tunately, because of the limited resolution of our calcu- 
lations we are not able to quantify this effect. The 
calculations illustrate the important point that magnetic 
buoyancy differs fundamentally from thermal convection, in 
which the effect of radiation is to cool the rising elements 
and reduce the convective flux. 

3.3. Saturation Phase 

At the end of the buoyant phase the disk undergoes a 
rather rapid transition to a state in which the initial, 
ordered field structure has evolved into a chaotic ensemble 
of patches (Fig. 6). However, while the patches are well 
defined and separated by field-free regions in the weak-seed 
case, in the strong-seed case they are diffuse, and the gas 
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Fig. 4.— Four different seeds in the inner disk are compared at an 
advanced stage of the buoyant phase (at the moment when the most 
evolved buoyant plume arrives at the surface of the disk), (a) Strong seed. 
(b-d) Weak seeds. Contours of constant density in the (r, z)-plane (solid 
curves) are shown with A log p = 0.28 and log p min = — 1 1.0 in all frames. 
Contours of constant ft with logarithmic spacing (broken curves) have log 
fi m in - 0.09, 0.45, 0.35, and 0.38 in ( a-d ), respectively, log /i max = 1.5 in (a), 
and log /J max - 1.8 in { b-d ). Arrows: Velocity vectors with F max = 1.30, 0.23, 
0.18, and 0.23 km s“ 1 in (a-d), respectively. Evolutionary times: 0.84, 4.15, 
4.32, and 2.22 yr in (a-d), respectively. The distance unit on the horizontal 
and vertical scales is 3.08 x 10 1 1 cm. 


between them is weakly magnetized. Further, in the strong- 
seed case the average saturation value of fi is close to unity, 
whereas in the weak-seed case it stays greater than five. In 
the strong-seed case the total magnetized area is larger, so 
that the difference in saturation values of the total field 
energy between the two cases approaches an order of mag- 
nitude. We conclude that at least up to evolutionary times 
indicated in Figure 6 the disk is able to preserve the 
memory of initial conditions. The explanation as to why the 
final total field energies are different is deferred to the end of 
§ 3.5. The disk states shown in Figure 6 are the final ones we 
were able to achieve. The simulations had to be stopped 
here because of uncontrolled field creation and annihilation 
on the grid-spacing scale. 

3.4. Energetics 

The monitoring of various forms of energy contained in 
the disk leads to further clarification of the physical pro- 



Fjg. 5. — Strong seeds evolving in the inner disk (a) with radiation 
transport and (b) adiabatically are compared at / = 0.82 yr. Curves, 
arrows, A log p, and log p mjn are the same as in Fig. 4. (a) and ( b ) log 
P min = —0.14 and —0.38, respectively; log (l mdX = 1.5 in both frames, (a) 
and (b) V mn = 1.28 and 0.53 km “ \ respectively. Distance unit : 3.08 x 10 l 1 
cm. 

cesses responsible for the evolutionary changes described in 
§§ 3. 1-3.3. As mentioned in § 3.1, in the buildup phase the 
azimuthal field grows at the expense of the orbital energy of 
the gas, which is driven by magnetic torques toward the 
center of the disk. At the same time, the poloidal field slowly 
decays as a result of numerical diffusion and/or reconnec- 
tion, and its energy E p decreases (Fig. 7). Thus, the total 
energy of the field £ f quickly becomes dominated by the 
energy of the azimuthal component, and, in accordance 
with equation (12), it grows roughly proportionally to t 2 . 
The accompanying growth of the kinetic energy K associ- 
ated with the motions in the meridional plane results 
mainly from the radial expansion of the seed loop. 

The total energy of the field keeps growing as ^ t 2 until 
the numerical loss of B ^ becomes significant. In the buildup 
phase, the only factor responsible for that loss is the auto- 



Fig. 6.— The final calculated states of (a) strong and (/>) weak seed 
evolving in the inner disk. Evolutionary times: (a) 4.90 yr and (7?) 1 1.90 yr. 
Curves, A log p, and log p mjn are the same as in Fig. 4. In frames (a) and (b) 
log p min = -0.09 and 0.74, log (i miX = 1.0 and 2.0, respectively. Distance 
unit: 3.08 x 10 1 1 cm. 
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Fig. 7.—Total magnetic energy (E,h total energy of the poloidal field 
(E p ), and total kinetic energy of motions in the meridional plane ( K ) as a 
function of time for the strong (subscript s) and weak (subscript w) seeds 
evolving in the inner disk. 


matic annihilation (reconnection) of oppositely oriented 
fluxes on the grid-spacing scale, which occurs whenever 
oppositely oriented patches of B $ field are advected into the 
same grid cell. In the early buildup phase, the interface 
dividing the seed loop into two halves with opposite signs of 
B# stays flat and smooth (Fig. la), there is little motion 
across it, and the annihilation rate is a small fraction of the 
generation rate. Only in the buoyant phase does the 
numerical reconnection become important enough to slow 
the growth of E r The primary cause of the increased 
numerical reconnection is vigorous buoyant flow, which 
causes the interface between the areas of oppositely oriented 
B $ to bend sharply, while the velocities perpendicular to it 
grow up to ~0.02c s in the weak seed case (Fig. 46) and 
~0.1c s in the strong seed case (Fig. 4a). The buoyant flow 
also causes the total area of the deformed loop to increase 
rapidly. As a result, the intensity of the azimuthal field 
decreases, owing to geometrical dilution of the azimuthal 
flux. This effect is probably more important than the 
numerical effect until the latest stages of the evolution. 

In the strong field case, the numerical loss and the geo- 
metrical dilution of acting together cause E t to drop by 
half an order of magnitude in less than ~0.5 yr between 
^0.75 and ~ 1.25 yr from the beginning of the simulation. 
The same effect is much gentler in the weak seed case, where 
E t decreases by a somewhat smaller amount between ~3.5 
and ^5.5 yr from the beginning of simulation (Fig. 7). In 
both cases, the poloidal energy E p reverses its downward 
trend and starts to increase during the epoch of declining 
E t , whereas the kinetic energy K markedly accelerates its 
growth before reaching a maximum value. As we already 
stressed in §§ 2.1 and 2.4, the null divergence of the poloidal 
field is always satisfied in our calculations (i.e., the total 
poloidal flux stays constant), which means that the 
observed increase of the poloidal field energy does not 
violate the axisymmetric antidynamo theorem (see, e.g., 
Parker 1979, § 18). E p grows simply because the poloidal 


field lines are wound more tightly; following the pictur- 
esque terminology of Hawley & Balbus (1992), we may say 
that the magnetized gas is kneaded by buoyancy-induced 
flows, with progressively smaller structures appearing in it 
as time goes on. It is worth noticing that both K and E p 
increase at the expense of the azimuthal field energy, which 
is first transformed into the kinetic energy of buoyant 
motions and then into the poloidal field energy as a result of 
the kneading phenomenon. 

At the moment the first plume arrives at the surface of the 
disk another sink of magnetic energy opens, originating 
from our condition of constant ambient pressure (§ 2.2). 
Any azimuthal field advected into the ambient medium by 
the residual mass transfer from the disk into the ambient 
medium is assumed to annihilate on the spot, and it is 
simply removed from the grid. Any poloidal field is prevent- 
ed from extending into the ambient medium through a pro- 
cedure in which the electromotive force (hereafter emf) in 
the induction equation is gradually reduced from its exact 
value just below the surface to zero in the ambient medium 
(note that, in the constrained transport method, the modifi- 
cations of emf do not violate the poloidal flux conservation 
condition). On the average, the surface loss of the azimuthal 
field is small (less than 1%) compared to the loss from 
numerical reconnection in the disk. Occasionally, however, 
it can exceed the latter because of a burstlike escape of B ^ 
into the ambient medium. The influence of emf reduction on 
the poloidal field evolution is also small as far as the total 
energy of the poloidal field is concerned; we checked that 
the E p curves obtained with and without reduction of the 
emf never differ by more than ^10%. However, this pro- 
cedure has a significant stabilizing effect on the evolution: if 
the emf is not reduced, the vigorous kneading just below the 
surface of the disk quickly brings the code to its limits. As 
far as the numerical reconnection is concerned, the extent to 
which it mimics real physical reconnection is difficult to 
estimate, since the true reconnection rate depends on the 
detailed geometry of the field, which is not known. This 
question is important, and it will be considered in the 
future. 

The epoch of E t decline comes to an end when a near 
equilibrium between field generation and destruction is 
reached in the saturation phase. B# is now generated from 
shearing of a highly nonuniform poloidal field which 
spreads over a larger and larger area of the disk. There is no 
well-defined interface separating regions with oppositely 
oriented B the motions are chaotic, and azimuthal field 
patches with different signs annihilate “on the spot,” at 
almost the same rate at which they are created, i.e., on a 
timescale of several orbital periods. It may be noted that 
this timescale is longer than the shortest possible timescale 
for physical reconnection of a field composed of thin 
patches, the latter being comparable to the orbital period 
(Rozyczka, Turner, & Bodenheimer 1995). As a result, the 
average intensity of the azimuthal field remains nearly con- 
stant, and the observed slow increase of E t is almost entirely 
caused by the increase of the total magnetized area. With 
the characteristic scale of magnetic field structure 
approaching the grid resolution, the numerical reconnec- 
tion and/or diffusion of the poloidal field also become more 
and more effective. As a result, in the saturation phase E p 
decreases slowly, with occasional short-term bursts owing 
to locally enhanced kneading. The fact that the kinetic 
energy of meridional motions also declines in the saturation 
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phase is a bit more difficult to explain. In order to elucidate 
this issue, let us recall that the increase of K at the end of the 
buoyant phase was mainly caused by energy injection on 
large scales (comparable to the characteristic size of the first 
strong buoyant plumes). Because the magnetic field has 
become more homogeneous, that process is much less effec- 
tive now (i.e., there is much less buoyant stirring on large 
scales), while on small scales K is partly converted into E p 
and partly dissipated by viscosity (see § 2.2) at about the 
same rate as before. 


3.5. Excitation Mechanism of the Buoyant Instability 

The initial conditions for our simulations are too compli- 
cated for a detailed linear stability analysis. However, even 
an approximate analysis based on a “submerged sharp 
boundary” case discussed by Parker (1979, p. 316) may be 
useful. To apply it, we will have to adopt two simplifying 
assumptions. First, we follow Parker by adopting the 
Boussinesq approximation, in which the effect of compress- 
ibility is neglected for the perturbed quantities. Second, we 
assume that the upper boundary of the magnetized region is 
originally entirely flat and infinite in extent (i.e., that the 
curvature of the initial poloidal field lines is unimportant), 
and that the generated azimuthal component is equivalent 
to a uniform field with straight lines perpendicular to the 
initial poloidal field lines. We shall also assume that the 
gravitational acceleration and total (gas plus magnetic) 
pressure are uniform in the vicinity of the upper boundary 
of the magnetized region. The stable lower boundary of the 
magnetized region will be entirely excluded from our con- 
siderations. 

To facilitate referencing to Parker’s work, let us adopt the 
same Cartesian coordinate system (x, y , z), and let us use the 
same symbols. Let the y-axis point in our azimuthal direc- 
tion, while x- and z-axes will point in our radial and vertical 
directions, respectively. Further, let the gravitational accel- 
eration g be oriented toward the midplane of the disk (i.e., in 
the negative z-direction). For the present analysis, we locate 
the interface between magnetized and field-free regions at 
z — 0. The magnetic field is perpendicular to g and confined 
to the half space z < 0, hereafter referred to as domain 1 (the 
field-free region is referred to as domain 2). Although the 
total pressure is continuous across the interface, the density 
in domain 2 is larger than that in domain 1 by an amount 
A p because of the drop in magnetic pressure. Neglecting 
changes in the sound speed c s , the total pressure 

p t = pci + (B 2 X + By)/Sn (14) 

in domain 1 maintains a balance with that in domain 2 
when 

A p = (B 2 X + B 2 )/Snc 2 , (15) 

where B x (our B r ) is approximately constant in time and B y 
(our Bj,) grows in time according to equation (12). 

Our first assumption results in 

V •» = (), (16) 

where v is the perturbed velocity. The rotational part of v 
only contributes to oscillatory Alfven waves and therefore is 
irrelevant here. In domain 1, the irrotational part of v may 
be expressed in terms of a potential function 

i? = — V4* , (17) 


with 

T = C(f) exp (ixk x + iyk y + Zyfkl + k 2 ) , (18) 

where C(r) is a function of f, k x , k y9 and k z are components of 
the wavevector Ac, and i is the imaginary unit. From the 
induction equation (4) we find the perturbed field com- 
ponents 

b x = (B,k x k, + B x klYTV. 
b y = (B y kj + B x k x k,YTV , (19) 

b z = (- iB, - iB x kjk 2 x + kJjW , (20) 

where V = f* 0 C(t')dt\ For domain 2, the potential function 
for the perturbed velocity is 

— S(t) exp (ixk x + iyk y — zj k 2 x -f k 2 ) . (21) 

The sign for the z dependence of the potential functions is 
chosen such that T and O vanish as z approaches - oo and 
+ oo, respectively. The location of the perturbed interface is 
given by a Lagrangian variable £(x, y, t) — \v z dt. The 
requirement for v x to be continuous across the interface 
implies that 

C(t) + S(t) = 0 (22) 

at z = Remembering that in domains 1 and 2 the pressure 
at the perturbed interface is lower than that at z = 0 by gp£ 
and g{p + Ap)c^, respectively, we can obtain from equation 
(3) the following condition for pressure balance at the inter- 
face (see also Parker 1 979, p. 3 1 8) : 

B B 

(bp 2 ) 0 = (bp 1)0 + ^ (b y ) 0 + ^ (b x ) 0 T g A pc, + A Q zz , 

(23) 

where dp is the pressure perturbation, subscripts 1 and 2 
refer to perturbations in magnetized and nonmagnetized 
domains, respectively, and subscript 0 denotes a value cal- 
culated at z = 0. The last term in equation (23) accounts for 
the nonvanishing difference in the zz component of the 
viscous stress tensor across the interface : 

= (2 p + Ap)v(k 2 x + /c 2 )T , (24) 


where v = rj/p and rj is given by equation (5). After appro- 
priate substitutions, equation (23) yields the amplitude 
equation 


(2 p + Ap)|^ + v(k 2 x + fc 2 )vj 

-[g^V^ JBrkx+ J’^ 


'Vdt' = 0 


(25) 


Since our numerical models are axisymmetric, in the follow- 
ing, let us focus on modes with k y = 0 which do not bend 
the lines of B y (i.e., our azimuthal field). We shall refer to 
them simply as interchange modes , keeping in mind, 
however, that only the azimuthal field lines are subject 
to pure interchange modes, while poloidal field lines may 
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still undulate. Under the substitutions x = Hr, k = Hk x , 
& = Cl 2 H, and c s = c s QH (where H is the disk half-thickness 
and CIH is the midplane value of the sound speed as given 
by the standard theory of thin accretion disks), equation 
(25) reduces to 


dC(x) ( 2c 2 ak 2 C(x) T 2 k 2 - k( 1 + 9t 2 /4)~ 

dx 3 + L 2/J 0 + 1 + 9t 2 /4 _ 


C(x r )(h f = 0 , 
(26) 


where /? 0 is the initial value of /?. 

The linear integrodifferential equation (26) for the ampli- 
tude C was solved numerically as an initial value problem 
with the help of a fourth-order Runge-Kutta method. To 
account for the displacement of the interface (i.e., the upper 
surface of the loop) from the midplane of the disk cl was set 
equal to 0.5. The results of integrations for a = 0.05 and 
a = 0.01 are plotted in Figure 8 as a function of k and t , 
with C normalized to unity at r = 0 for all k. In each panel 
the curves are separated by equal intervals of time (At = 3 
and 10 for strong and weak seed models, respectively). The 
highlighted curves in the a = 0.05 panels correspond to the 
beginning of the buoyant phase as determined from our 
simulations (epoch t = 15 for the strong seed and t = 80 for 
the weak seed). Note that, at these epochs, maximum C is 
attained for k ~~ 5-6, which agrees reasonably well with the 
wavelengths actually excited (Figs. 3 and 4). 


Figure 8 also indicates that short-wavelength modes 
grow at a significantly lower rate than the most unstable 
mode. This finding runs contrary to the original results of 
Parker, according to which, in the submerged sharp bound- 
ary case, shorter wavelengths have consistently higher 
growth rates. This difference is a result of the presence of 
both the poloidal field and the viscosity (factors which were 
absent from Parker’s analysis). In equation (26), the viscous 
(i.e., the second) term clearly damps the growth of perturbed 
quantities. The effectiveness of the viscous damping may be 
estimated with the help of the dimensionless damping time- 
scale 


*d = 



(27) 


In our numerical simulations, z d approaches unity (i.e., the 
dynamical timescale) for wavelengths comparable to //, 
while perturbations with X < H are damped on timescales 
shorter than the dynamical timescale. A similar stabilizing 
role is played by the tension of the B x (our B r ) lines. In order 
to elucidate this effect, let us set v = 0 and define a charac- 
teristic timescale of the amplitude growth 


C 

9 dC/dz 


(28) 


Despite the fact that z g decreases with t initially and 






Fig. 8. Results from linear stability analysis. The amplitude of the perturbation is plotted as a function of wavenumber k — 2 nH/X at several evolution- 
ary times, where H is the half-thickness of the disk. Each frame is labeled with the assumed initial ratio (fi) of gas pressure to magnetic pressure and with the 
viscosity parameter (a). The times for the cases fi = 3 x 10 2 start at Qf = 3 (lowest curves ) and increase upward in increments of Ut — 3. The times for the 
cases p = 3 x 10 4 start at Of = 10 (lowest curves) and increase upward in increments of Qt = 10. Heavy lines : time of onset of instability in the numerical 
simulations. 
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approaches a constant for large t, let us approximate 

“ do • 

Based on this approximation, equation (26) reduces to 
k( 1 + 9t 2 /4) 


1 Ik 2 
— + 

*0 


20 o + 1 + 9 t 2 /4 Tg ° ’ 


(29) 


(30) 


yielding the following stability condition in the absence of 
viscosity : 




(31) 


A comparison of the corresponding left and right panels in 
Figure 8 shows that for the parameters considered here, the 
viscous effects on the stability of the buoyant modes may 
dominate those resulting from the tension of B r lines. This 
result, however, is not absolutely certain. Also, the depen- 
dence of the most unstable wavelength on the initial seed 
strength, clearly visible in the simulations, is detected only 
marginally by the linear stability analysis (compare the two 
left panels with a = 0.05 in Fig. 8). These two facts indicate 
that additional calculations and a more detailed analysis 
would be necessary to clearly assess the role of both factors. 
However, such an analysis is unlikely to change the main 
conclusion of the present section: the plumes observed in 
our simulations originate from interchange modes acting 
predominantly on the azimuthal field, and mediated by the 
presence of the poloidal field and nonnegligible viscosity. 
Moreover, in accordance with the simulations, the present 
analysis predicts that the most unstable wavelength should 
scale with the half-thickness of the disk. Based on this 
analysis, we can also provide an explanation for the differ- 
ence in the final total field energies mentioned in § 3.3. Mea- 
suring the logarithmic derivative of the perturbation 
amplitude d In C/dz, we find that it first approaches unity at 
the peak of the highlighted curves in the left panels of 
Figure 8. The corresponding timescale for amplitude 
growth z g also approaches unity, which means that the per- 
turbations begin to grow on a dynamical timescale, and a 
transition from the buildup to the buoyant phase is initial- 
ized. At the moment when z g drops to unity, (i predicted by 
the linear analysis approaches ~ 0.6 and ^2.1 for the 
strong and weak seed, respectively. Shortly afterward, 
strong buoyant motions begin to stir the disk, and the 
reconnection rate goes up, preventing further growth of B 
Thus, the field saturates at a strength not much different 
from the one it had when z g approached unity for the fastest 
growing mode. 


4. SUMMARY AND DISCUSSION 

The simulations reported in the preceding sections fol- 
lowed the evolution of a seed poloidal field in an axisym- 
metric accretion disk. We found that the azimuthal 
component of the field B $ grows as a result of shear until the 
thermal equipartition value of (i = 1 is approached. At the 
same time, the seed spreads radially (toward both the center 
of the disk and its outer edge) as a result of the magneto- 
rotational instability (hereafter MRI). Within the seed, the 
angular momentum of the disk gas is removed by magnetic 
torques from mass elements at smaller orbital radii, and it is 
added to mass elements at larger orbital radii. The growth 
of B ^ is stopped by vigorous buoyant motions leading even- 


tually to partial ejection of the field from the disk, fragmen- 
tation of the seed, and the establishment of a turbulent 
velocity field. The linear stability analysis indicates that the 
buoyancy effects are a result of interchange modes of the 
Parker instability (Parker 1979, § 13) acting on the azi- 
muthal component of the field, with the effectively excited 
wavelength depending on the poloidal seed strength and the 
effective viscosity assumed for the disk gas. At later evolu- 
tionary times, a quasi-stationary state is reached, in which 
the total magnetic energy of the disk slowly grows because 
of radial expansion of the magnetized region, while the 
average values of field components either remain roughly 
constant or decrease slowly owing to numerical reconnec- 
tion. In that state, the poloidal field assumes a patchy struc- 
ture, with the medium between patches being either 
field-free or more weakly magnetized than the patches 
themselves. (Let us note that a patch of the poloidal field, 
for example, one of those plotted in Fig. 6 b, is also a merid- 
ional cross section of a bundle of toroidal field lines). 

In our models, both the field intensity at which buoyancy 
begins to effectively operate and the saturation intensity 
decrease with decreasing initial seed field strength. The 
dependence, however, is rather weak (see § 3.3). One may 
safely assume that in the presence of B ^ -+ B r feedback the 
seed fields would be amplified, and it is likely that a satura- 
tion intensity independent of the initial seed strength (i.e., 
defined solely by the disk’s microphysics) would be 
achieved. Unless the three-dimensional effects or physical 
dissipation of magnetic fields prove to be very important, 
on the basis of our simulations, it is rather difficult to 
imagine that the field would saturate at P P 1. Thus, at least 
in the realm of ideal MHD, we expect the buoyant dynamo 
to generate rather strong fields. On the other hand, our 
calculations indicate that P — 1 is an absolute limit to field 
amplification, since already at % 1 buoyant motions are 
so vigorous that the field begins to escape from the disk. 
The same limiting field strength was obtained by Galeyev et 
al. (1979), who discussed a classical turbulent disk dynamo. 
Let us note that at P ^ 1 the Alfven velocity, which may be 
regarded as a reasonable approximation to the buoyancy- 
generated velocities (Meyer & Meyer-Hofmeister 1983; 
Tout & Pringle 1995), is comparable to the velocity of 
sound c s . The characteristic timescale for the field escape is 
then only H/c s ^ 1/Q. 

From our numerical results we can now estimate the 
effective viscosity parameter a introduced by Shakura & 
Sunyaev (1973). For simplicity, we shall assume that in the 
saturation state the field is indeed collected in rather well- 
defined patches (see Fig. 6b), and that each patch evolves in 
the same way as the initial seed loop, experiencing its own 
buildup phase (however, with a value of B r which does not 
necessarily have to be the same as that of the initial seed). 
We further assume that the buildup phase always ends at 
p = 1, neglecting the weak dependence of saturation inten- 
sity on B r . In other words, we assume that the field saturates 
when the magnetic pressure generated by the azimuthal 
component of the field becomes equal to the thermal pres- 
sure. If the initial j? of a typical patch is significantly smaller 
than unity, and if B r stays roughly constant throughout the 
buildup phase, then it follows from equation (12) that the 
time t h needed to make the azimuthal field effectively 
buoyant, i.e., to achieve the saturation value of 

B 4, = (W 8n P ’ 


(32) 
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is given by 


h 3 QB r ' J ^ P 


(33) 


During that time the accumulated magnetic torque per unit 
mass results in a specific angular momentum change of 



(34) 


resulting in an average magnetic torque 


where 


AJ 2 cl 


(35) 


Pr = 


Snpc; 

B] 


Comparing j m with the viscous “ a-torque,” 


J, = j ac * 


(36) 


(37) 


we get 


a 


4 1 

3 TiT 


(38) 


Setting P r = p 0 in order to apply this general formula to our 
original strong and weak seed fields, we obtain an effective a 
of 0.06 and 0.006, respectively. We see that in the strong 
field case a (coincidental) agreement is found between the 
assumed and the effectively generated value of oc. In two- 
dimensional calculations, the same a-values are applicable 
during the saturation stage, because B r approximately pre- 
serves its initial value. However, in three-dimensional calcu- 
lations B r would be expected to increase, resulting in an 
increase in the effective value of a during the evolution. 

The enforced axial symmetry of our simulations prevents 
the transformation of into the poloidal components B z 
and B r . However, in three dimensions such a transform- 
ation is unavoidable, since the Coriolis force acting on 
buoyant elements results in helical twisting of the azimuthal 
field lines. The feedback between B ^ and B r is a necessary 
requirement for the dynamo operation. Given the robust- 
ness of the buoyant motions illustrated by our simulations, 
it seems that such a dynamo, powered solely by shear and 
hydromagnetic instabilities, should be a vital component of 
the physical processes governing the structure and evolu- 
tion of accretion disks. Moreover, its operational range 
would extend into convectively stable regions of the disks, 
where field generation mechanisms that rely on an under- 
lying turbulence cannot be applied. In a convectively stable 
region, the amplification of the seed field would not be ham- 
pered by the nonlinear effects discussed by Kulsrud & 
Anderson (1992) and Vainshtein & Cattaneo (1992), since 
the turbulence would appear only after a strong azimuthal 
field has been generated. 

In an idealized scenario, the dynamo would be activated 
whenever a radial seed field appeared in the disk. The 
primary energy (the orbital energy of the disk gas) would be 
transformed into the azimuthal field energy, making the 
seed increasingly buoyant. The buoyancy effect in turn 
would transform the azimuthal field energy into kinetic 


energy of radial and vertical motions of the disk gas. These 
motions alone, without the help of the Coriolis effect, would 
be able to turn their own energy partly into the energy of 
the poloidal field through the kneading phenomenon dis- 
cussed in § 3.4. At the same time, the vertical motions com- 
bined with the Coriolis effect would regenerate and/or 
amplify the initial seed. Dissipation of kinetic and magnetic 
energies within the disk would generate heat, resulting in an 
enhancement of buoyant motions by radiative energy 
inflow into the rising elements (see § 3.2 and Stella & 
Rosner 1984). Because of the MRI, the whole dynamo 
process would be necessarily accompanied by angular 
momentum transfer from lower to higher orbits and by an 
accretion flow. Thus, the expansion of an initially localized 
dynamo process would ensue, which could only be stopped 
in the regions where coupling between the field and the disk 
gas is inefficient owing to a very low degree of ionization. 

However, many questions have to be answered before the 
above-outlined scenario could be adopted. Probably the 
most important one among them concerns the long-term 
evolution of the buoyant dynamo over periods comparable 
to, or longer than, the accretion timescale. One may suspect 
that the same buoyancy that is responsible for magne- 
tization of the disk will slowly remove the field from the 
vicinity of the equatorial plane and drive the magnetized 
area toward the surface, eventually leading to complete 
demagnetization of the disk gas. A counterargument is that 
buoyancy does not operate at z = 0, and thus the removal 
of the magnetic field from the equatorial plane would not 
occur. Extensive numerical research is required to deter- 
mine whether a strong equilibrium field can be maintained 
in the long term. Further, three-dimensional effects may not 
only support the dynamo action but also counteract it. For 
example, BGH find that weak azimuthal fields are also 
prone to the MRI. Thus, a suspicion arises that the buildup 
phase could be substantially shortened, saturating the 
dynamo at a much lower field intensity than suggested 
above. On the other hand, BGH point out that as time 
advances the wavelength of the dominant azimuthal mode 
of the MRI becomes longer. This would mean that the field 
is able to organize itself into more coherent structures 
which may eventually evolve into nearly axisymmetric 
patches. Another three-dimensional effect whose influence 
on the efficiency of the buoyant dynamo remains unclear is 
the undulating mode of the Parker instability, as a result of 
which the azimuthal field lines bend vertically. Again, exten- 
sive numerical research is needed to elucidate these issues. 
Finally, one may ask how the buoyant dynamo would 
operate in the presence of a turbulent velocity field not 
related to magnetic buoyancy (e.g., excited by thermal 
convection), and how it would react to an externally con- 
nected field threading the disk. While the answer to the 
latter qestion may be connected to the BGH results that 
indicate enhanced magnetoturbulence caused by the MRI, 
the interaction between convective flows and magnetic 
buoyancy effects seems to be too difficult a problem to be 
tackled at the present state of research. 
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